text:
# MIT License (MIT)
#
# Copyright (c) 2015 Christoph Feinauer
#
# 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.

# For questions and suggestions please use the Github page https://github.com/christophfeinauer/PdbTool

module PdbTool
using Compat

	macro spath()
	        return dirname(Base.source_path()) * "/"
	end

	KRED="\x1B[31m"
	KGRN="\x1B[32m"
	KRES="\033[0m"

	#include(@spath()*"../test/"*"runtests.jl")
	include("EXT_TEST.jl")
	include("aminoAcidDict.jl")
	using .EXT_TEST

	######################################################################
	# TYPE DEFINITIONS
	######################################################################
	type Atom
	        identifier::Int64
		@compat coordinates::Tuple{Float64,Float64,Float64}
	end
	type Residue
		aminoAcid::AbstractString
	        atom::Dict{AbstractString,Atom}
		pdbPos::Int64
		alignmentPos::Int64
		identifier::AbstractString
		naccess::Float64
		naccess_rel::Float64
		Residue(aa,pp,id)=new(aa,Dict{AbstractString,Atom}(),pp,-1,id,-1.0,-1.0)
	end
	type Strand
		identifier::Int64
		startRes::AbstractString
		endRes::AbstractString
		sense::Int64
		bondThis::AbstractString
		bondPrev::AbstractString
		Strand(identifier::Int64,startRes::AbstractString,endRes::AbstractString,sense::Int64,bondThis::AbstractString,bondPrev::AbstractString)=new(identifier::Int64,startRes::AbstractString,endRes::AbstractString,sense::Int64,bondThis::AbstractString,bondPrev::AbstractString)
	end
	type Sheet
		identifier::AbstractString
		strand::Dict{Int64,Strand}
		Sheet(identifier::AbstractString)=new(identifier,Dict{Int64,Strand}())
	end
	type Helix
		startRes::AbstractString
		endRes::AbstractString
		identifier::AbstractString
		Helix(sR,eR,id)=new(sR,eR,id)
	end
	type Chain
	        residue::Dict{AbstractString,Residue}
	        length::Int64
		mappedTo::AbstractString
		identifier::AbstractString
		align::Dict{Int64,Residue}
		helix::Dict{AbstractString,Helix}
		sheet::Dict{AbstractString,Sheet}
		isRNA::Bool
	        Chain()=new(Dict{AbstractString,Residue}(),0,"","",Dict{Int64,Residue}(),Dict{AbstractString,Helix}(),Dict{AbstractString,Sheet}(),false)
	end
	type Pdb
	        chain::Dict{AbstractString,Chain}
	        pdbName::AbstractString
	        fileName::AbstractString
	        Pdb()=new(Dict{AbstractString,Chain}(),"","")
	end
	include("show_methods.jl")

	######################################################################
	# FUNCTION:		 parsePdb
	######################################################################
	function parsePdb(pdbFile::AbstractString)
		!isfile(pdbFile) && error("File not found")
		pdb=Pdb()
		for l in eachline(open(pdbFile))
			if l[1:4]=="ATOM" && (l[17]==' ' || l[17]=='A')
				ch=strip(l[22:22])
				res=strip(l[23:27])
				if !haskey(pdb.chain,ch)
					pdb.chain[ch]=Chain()
					pdb.chain[ch].identifier=ch
					if length(lstrip(l[18:20]))<3
						pdb.chain[ch].isRNA=true
						println("Chain $ch treated as RNA")
					end
				end
				if !haskey(pdb.chain[ch].residue,res)
					pdb.chain[ch].length += 1
					pdb.chain[ch].residue[res]=Residue(l[18:20],pdb.chain[ch].length,res)
				end
				pdb.chain[ch].residue[res].atom[strip(l[13:16])]=Atom(parse(Int,l[7:11]),(parse(Float64,l[31:38]),parse(Float64,l[39:46]),parse(Float64,l[47:54])))
				if pdb.chain[ch].isRNA
					for res in values(pdb.chain[ch].residue)
						res.aminoAcid=lstrip(res.aminoAcid)
					end
				end

			end
			if l[1:5]=="HELIX"
				ch=strip(l[20:20])
				hel=strip(l[12:14])
				if !haskey(pdb.chain,ch)
					pdb.chain[ch]=Chain()
					pdb.chain[ch].identifier=ch
					if length(lstrip(l[18:20]))<3
						pdb.chain[ch].isRNA=true
						println("Chain $ch treated as RNA")
					end
				end
				if haskey(pdb.chain[ch].helix,hel)
					error("Found the same helix twice and panicked")
				end
				pdb.chain[ch].helix[hel]=Helix(strip(l[22:26]),strip(l[34:38]),hel)
			end
			if l[1:5]=="SHEET"
				if l[22:22]!=l[33:33]
					warn("PDB contains sheets across several chains - ignoring them for now!")
					continue
				end
				ch=strip(l[22:22])
				sh=strip(l[12:14])
				str=parse(Int64,strip(l[8:10]))
				if !haskey(pdb.chain,ch)
					pdb.chain[ch]=Chain()
					pdb.chain[ch].identifier=ch
					if length(lstrip(l[18:20]))<3
						pdb.chain[ch].isRNA
						println("Chain $ch treated as RNA")
					end
				end
				if !haskey(pdb.chain[ch].sheet,sh)
					pdb.chain[ch].sheet[sh]=Sheet(sh)
				end
				if haskey(pdb.chain[ch].sheet[sh].strand,str)
					error("Found the same strand twice and panicked")
				end
				sense=parse(Int64,l[39:40])
				if sense!=0
					pdb.chain[ch].sheet[sh].strand[str]=Strand(str,strip(l[23:27]),strip(l[34:38]),sense,strip(l[51:55]),strip(l[66:70]))
				else
					pdb.chain[ch].sheet[sh].strand[str]=Strand(str,strip(l[23:27]),strip(l[34:38]),sense,"","")
				end
			end

		end
		return pdb;
	end

	######################################################################
	# FUNCTION:		 atomDist
	######################################################################
	function atomDist(atom1::Atom, atom2::Atom)
			return sqrt((atom1.coordinates[1]-atom2.coordinates[1])^2 + (atom1.coordinates[2] - atom2.coordinates[2])^2 + (atom1.coordinates[3] - atom2.coordinates[3])^2)
	end

	######################################################################
	# FUNCTION:		 residueDist
	######################################################################
	function residueDist(res1::Residue, res2::Residue; distType="heavyMin")
		if distType=="heavyMin"
			d=Inf
			for atom1 in values(res1.atom)
				for atom2 in values(res2.atom)
					d=min(d,atomDist(atom1,atom2))
				end
			end
			return d
		end
		if distType=="ca"
			if haskey(res1.atom,"CA") && haskey(res2.atom,"CA")
				return atomDist(res1.atom["CA"],res2.atom["CA"])
			else
			 	error("error calculating distances: Residues do not have CA entries")
			end
		end
		error("Error calculating distances: supported types are heavyMin or ca")
	end

	######################################################################
	# FUNCTION:		 chainDist
	######################################################################
	function chainDist(ch1::Chain, ch2::Chain,out="undef")
		if out=="undef"
			for r1=sort([k for k in keys(ch1.residue)])
				for r2=sort([k for k in keys(ch2.residue)])
					@printf("%s %s %f\n", r1, r2, residueDist(ch1.residue[r1],ch2.residue[r2]))
				end
			end
		else
			fid=open(out,"w")
			for r1=sort([k for k in keys(ch1.residue)])
				for r2=sort([k for k in keys(ch2.residue)])
					@printf(fid,"%s %s %f\n", r1, r2, residueDist(ch1.residue[r1],ch2.residue[r2]))
				end
			end
		end
	end

	######################################################################
	# FUNCTION:		 interAlignDist
	######################################################################
	function interAlignDist(ch1::Chain, ch2::Chain;out="undef")
		if ch1.mappedTo==""
			error("chain 1 not mapped")
		end
		if ch2.mappedTo==""
			error("chain 2 not mapped")
		end
		LENG1=getHmmLength(ch1.mappedTo)
		LENG2=getHmmLength(ch2.mappedTo)
		completeAlign=Dict{Int64,Residue}()
		for i in keys(ch1.align)
			completeAlign[i]=ch1.align[i]
		end
		for i in keys(ch2.align)
			completeAlign[i+LENG1]=ch2.align[i]
		end
		ind=sort([x for x in keys(completeAlign)])
		if out=="undef"
			for i=1:length(ind)
				for j=(i+1):length(ind)
					@printf("%s %s %f\n", ind[i], ind[j], residueDist(completeAlign[ind[i]],completeAlign[ind[j]]))
				end
			end
		else
			fid=open(out,"w")
			for i=1:length(ind)
				for j=(i+1):length(ind)
					@printf(fid,"%s %s %f\n", ind[i], ind[j], residueDist(completeAlign[ind[i]],completeAlign[ind[j]]))
				end
			end
			close(fid)
		end
	end

	######################################################################
	# FUNCTION:		 chainSeq
	######################################################################
	function chainSeq(chain::Chain)
		# Unordered pdbSeq
		if !chain.isRNA
			pdbSeq=[aminoAcidDict[chain.residue[k].aminoAcid] for k in keys(chain.residue)]
		else
			pdbSeq=[chain.residue[k].aminoAcid for k in keys(chain.residue)]
		end

		ind=sortperm([chain.residue[k].pdbPos for k in keys(chain.residue)])
		return join(pdbSeq[ind])
	end

	######################################################################
	# FUNCTION:		 mapChainToHmm
	######################################################################
	function mapChainToHmm(chain::Chain,hmmFile::AbstractString)
		# Check if the chain already has a mapping - and delete it if yes
		if chain.mappedTo!=""
			println("chain $(chain.identifier) already has a mapping.")
		end

		# Actual mapping
		pdbSeq=chainSeq(chain)
		tempFile=tempname()
		fid=open(tempFile,"w")
		@printf(fid,">temp\n%s",pdbSeq)
		close(fid)
		if !chain.isRNA
			if !EXT_TEST_hmmalign()
				error("cannot run hmmalign - please check that it is on the path")
			end
                    @compat run(pipeline(`hmmalign $hmmFile $tempFile`,"$tempFile.out"))
		else
			if !EXT_TEST_cmsearch()
				error("cannot run cmsearch - please check that it is on the path")
			end

                    @compat run(pipeline(`cmsearch -A $tempFile.out $hmmFile $tempFile`,DevNull))
		end

		st2fa("$tempFile.out";oFile=tempFile)
	        @compat align=collect(split(readall(tempFile),'\n'))

		rm("$tempFile")
		rm("$tempFile.out")
		pdbIndices=find([align[2][x]!='-' for x=1:length(align[2])])
		cleanIndices=find(![islower(align[2][x]) for x=1:length(align[2])])

		fakeAlign2pdb=-ones(Int64,length(align[2]))
		@compat fakeAlign2pdb[pdbIndices]=collect(1:length(pdbSeq))
		align2pdb=fakeAlign2pdb[cleanIndices]

		for k in keys(chain.residue)
			if chain.residue[k].pdbPos > 0
				x=find(align2pdb.==chain.residue[k].pdbPos)
				if length(x)>1
					error("found several pdb positions")
				end
				if length(x)==1
					chain.residue[k].alignmentPos=x[1]
					chain.align[x[1]]=chain.residue[k]
				end
			end
		end
		chain.mappedTo=hmmFile
	end

	######################################################################
	# FUNCTION:		 alignSeqToHmm
	######################################################################
	function alignSeqToHmm(seq::AbstractString,hmmFile::AbstractString)

		# Actual mapping
		tempFile=tempname()
		fid=open(tempFile,"w")
		@printf(fid,">temp\n%s",seq)
		close(fid)
		if !EXT_TEST_hmmalign()
			error("cannot run hmmalign - please check that it is on the path")
		end
                @compat run(pipeline(`hmmalign $hmmFile $tempFile`,"$tempFile.out"))

		st2fa("$tempFile.out";oFile=tempFile)
	        @compat align=collect(split(readall(tempFile),'\n'))

		rm("$tempFile")
		rm("$tempFile.out")
		return filter(x->!islower(x),align[2])
	end

	######################################################################
	# FUNCTION:		 mapSeqToHmm
	######################################################################
	function mapSeqToHmm(seq::AbstractString,hmmFile::AbstractString)

		# Actual mapping
		tempFile=tempname()
		fid=open(tempFile,"w")
		@printf(fid,">temp\n%s",seq)
		close(fid)
		if !EXT_TEST_hmmalign()
			error("cannot run hmmalign - please check that it is on the path")
		end
                @compat run(pipeline(`hmmalign $hmmFile $tempFile`,"$tempFile.out"))

		st2fa("$tempFile.out";oFile=tempFile)
	        @compat align=collect(split(readall(tempFile),'\n'))

		rm("$tempFile")
		rm("$tempFile.out")
		seqIndices=find([align[2][x]!='-' for x=1:length(align[2])])
		cleanIndices=find(![islower(align[2][x]) for x=1:length(align[2])])

		fakeAlign2seq=-ones(Int64,length(align[2]))
		@compat fakeAlign2seq[seqIndices]=collect(1:length(seq))
		align2seq=fakeAlign2seq[cleanIndices]

		hmm2seq = Dict{Int64,Int64}()
		for (hmm_ind,seq_ind) in enumerate(align2seq)
			seq_ind == -1 && continue
			hmm2seq[hmm_ind] = seq_ind
		end
		return hmm2seq
	end

	######################################################################
	# FUNCTION:		 mapChainToHmmLegacy
	######################################################################
	function mapChainToHmmLegacy(chain::Chain,hmmFile::AbstractString)
		# Check if the chain already has a mapping - and delete it if yes
		if chain.mappedTo!=""
			println("chain $(chain.identifier) already has a mapping.")
		end

		# Actual mapping
		pdbSeq=chainSeq(chain)
		tempFile=tempname()
		fid=open(tempFile,"w")
		@printf(fid,">temp\n%s",pdbSeq)
		close(fid)
		if !chain.isRNA
			if !EXT_TEST_hmmsearch()
				error("cannot run hmmsearch - please check that it is on the path")
			end
                    @compat run(pipeline(`hmmsearch -A $tempFile.out $hmmFile $tempFile`, DevNull))
		else
			if !EXT_TEST_cmsearch()
				error("cannot run cmsearch - please check that it is on the path")
			end
		    @compat run(pipeline(`cmsearch -A $tempFile.out $hmmFile $tempFile`, DevNull))

		end
		st2fa("$tempFile.out";oFile=tempFile)
		align=[split(readall(tempFile),'\n')]
		rm("$tempFile")
		rm("$tempFile.out")
		(pdbStart,pdbStop)=parse(Int64,matchall(r"\d+",align[1]))
		pdbIndices=find([align[2][x]!='-' for x=1:length(align[2])])
		cleanIndices=find(![islower(align[2][x]) for x=1:length(align[2])])
		fakeAlign2pdb=-ones(Int64,length(align[2]))
		fakeAlign2pdb[pdbIndices]=[pdbStart:pdbStop]
		align2pdb=fakeAlign2pdb[cleanIndices]
		for k in keys(chain.residue)
			if chain.residue[k].pdbPos > 0
				x=find(align2pdb.==chain.residue[k].pdbPos)
				if length(x)>1
					error("found several pdb positions")
				end
				if length(x)==1
					chain.residue[k].alignmentPos=x[1]
					chain.align[x[1]]=chain.residue[k]
				end
			end
		end
		chain.mappedTo=hmmFile
	end

	######################################################################
	# FUNCTION:		 intraAlignDist
	######################################################################
	function intraAlignDist(chain::Chain;out="distMat")
		if out=="distMat"
			if chain.mappedTo==""
				error("chain has no mapping")
			end
			LENG=getHmmLength(chain.mappedTo)
			distMat=-ones(LENG,LENG)
			for k1 in keys(chain.residue)
				for k2 in keys(chain.residue)
					if k1==k2
						continue
					end
					r1=chain.residue[k1]
					r2=chain.residue[k2]
					if r1.alignmentPos > 0 && r2.alignmentPos > 0
						distMat[r1.alignmentPos,r2.alignmentPos] = residueDist(r1,r2)
					end
				end
			end
		end
		return distMat
	end

	######################################################################
	# FUNCTION:		 makeIntraRoc
	######################################################################
	@compat function makeIntraRoc(score::Array{Tuple{Int64,Int64,Float64},1},chain::Chain;sz=200,cutoff::Float64=8.0,out::AbstractString="return",pymolMode::Bool=false,minSeparation::Int64=4)
		if chain.mappedTo==""
			error("chain has no mapping")
		end
		if out!="return"
			fid=open(out)
		end
		if !pymolMode
			roc=Array(Tuple{AbstractString,AbstractString,Float64,Float64},0)
		else
			roc=Array(Tuple{AbstractString,AbstractString,Int64,Float64},0)
		end
			s::Int64=0
			i::Int64=0
			hits::Int64=0
			while s<sz && i<length(score)
				i+=1
				if abs(score[i][1]-score[i][2]) <= minSeparation
					continue
				end
				if haskey(chain.align,score[i][1]) && haskey(chain.align,score[i][2])
					s+=1
					if residueDist(chain.align[score[i][1]],chain.align[score[i][2]])<cutoff
						hits+=1
						if pymolMode
							x=1
						else
							x=hits/s
						end
					else
						x=hits/s
						if pymolMode
							x=0
						end
					end
					push!(roc,(chain.align[score[i][1]].identifier,chain.align[score[i][2]].identifier,x,score[i][3]))
				end
			end
		return roc
	end

	######################################################################
	# FUNCTION:		 filterInterScore
	######################################################################
	@compat function filterInterScore(score::Array{Tuple{Int64,Int64,Float64},1},chain1::Chain,chain2::Chain;sz=200,cutoff::Float64=8.0,out::AbstractString="return")

		# Check if mapping is existent
		if chain1.mappedTo == ""
				error("chain 1 has no mapping")
		elseif chain2.mappedTo == ""
				error("chain 2 has no mapping")
		end
		LENG1=getHmmLength(chain1.mappedTo)
		LENG2=getHmmLength(chain2.mappedTo)

		if out=="return"
			newScore=Array((AbstractString,AbstractString,Float64),sz)
			s::Int64=0
			i::Int64=0
			hits::Int64=0
			positives::Int64=0
			while s<sz && i<size(score,1)
				i+=1
				if score[i][1] <= LENG1 && score[i][2] > LENG1
					ind1=score[i][1]
					ind2=score[i][2]-LENG1
				elseif score[i][2] <= LENG1 && score[i][1] > LENG1
					ind1=score[i][2]
					ind2=score[i][1]-LENG1
				else
					continue
				end
				if haskey(chain1.align,ind1) && haskey(chain2.align,ind2)
					s+=1
					id1=chain1.align[ind1].identifier
					id2=chain2.align[ind2].identifier
					newScore[s]=(id1,id2,score[i][3])
				end
			end
			return newScore
		end
	end
	@compat function filterInterScore(score::Array{Tuple{Int64,Int64,Float64},1},hmm1::AbstractString,hmm2::AbstractString;sz=200,cutoff::Float64=8.0,out::AbstractString="return")

		# Check if mapping is existent
		LENG1=getHmmLength(hmm1)
		LENG2=getHmmLength(hmm2)

		if out=="return"
			interScore=Array(Tuple{Int64,Int64,Float64},0)
			score1=Array(Tuple{Int64,Int64,Float64},0)
			score2=Array(Tuple{Int64,Int64,Float64},0)

			s::Int64=0
			i::Int64=0
			hits::Int64=0
			positives::Int64=0
			for i=1:length(score)

				if score[i][1] <= LENG1 && score[i][2] > LENG1
					ind1=score[i][1]
					ind2=score[i][2]-LENG1
					push!(interScore,(ind1,ind2,score[i][3]))

				elseif score[i][2] <= LENG1 && score[i][1] > LENG1
					ind1=score[i][2]
					ind2=score[i][1]-LENG1
					push!(interScore,(ind1,ind2,score[i][3]))

				elseif score[i][1] <= LENG1 && score[i][2] <= LENG1
					ind1=score[i][1]
					ind2=score[i][2]
					push!(score1,(ind1,ind2,score[i][3]))

				elseif score[i][1] > LENG1 && score[i][2] > LENG1
					ind1=score[i][1]-LENG1
					ind2=score[i][2]-LENG1
					push!(score2,(ind1,ind2,score[i][3]))
				end
			end
			return interScore,score1,score2
		end
	end

	######################################################################
	# FUNCTION:		 makeInterRoc
	######################################################################
	@compat function makeInterRoc(score::Array{Tuple{Int64,Int64,Float64},1},chain1::Chain,chain2::Chain;sz=200,cutoff::Float64=8.0,out::AbstractString="return",pymolMode::Bool=false,naccessRatio::Float64=1.0)

		# Check if mapping is existent
		if chain1.mappedTo == ""
				error("chain 1 has no mapping")
		elseif chain2.mappedTo == ""
				error("chain 2 has no mapping")
		end
		LENG1=getHmmLength(chain1.mappedTo)
		LENG2=getHmmLength(chain2.mappedTo)

		# Get naccess cutoff if necessary
		if naccessRatio<1.0
			naList1=zeros(LENG1);
			naList2=zeros(LENG2);
			i=1;
			for r1 in values(chain1.align)
				naList1[i]=r1.naccess;
				i+=1;
			end
			i=1;
			for r2 in values(chain2.align)
				naList2[i]=r2.naccess;
				i+=1;
			end
			naList1=sort(naList1,rev=true);
			na1Cutoff=naList1[parse(Int64,round(LENG1*naccessRatio))];
			println(na1Cutoff)
			naList2=sort(naList2,rev=true);
			na2Cutoff=naList2[parse(Int64,round(LENG2*naccessRatio))];
			println(na2Cutoff)
		end

		if out=="return"
			roc=Array(Tuple{AbstractString,AbstractString,Float64,Float64},0)
			s::Int64=0
			i::Int64=0
			hits::Int64=0
			positives::Int64=0
			while s<sz && i<size(score,1)
				i+=1
				if score[i][1] <= LENG1 && score[i][2] > LENG1
					ind1=score[i][1]
					ind2=score[i][2]-LENG1
				elseif score[i][2] <= LENG1 && score[i][1] > LENG1
					ind1=score[i][2]
					ind2=score[i][1]-LENG1
				else
					continue
				end
				if haskey(chain1.align,ind1) && haskey(chain2.align,ind2)
					if naccessRatio<1.0
						if chain1.align[ind1].naccess<na1Cutoff || chain2.align[ind2].naccess<na2Cutoff
							continue;
						end
					end
					s+=1
					id1=chain1.align[ind1].identifier
					id2=chain2.align[ind2].identifier
					if residueDist(chain1.align[ind1],chain2.align[ind2])<cutoff
						hits+=1
						if pymolMode
							hits=s
						end
						push!(roc,(id1,id2,hits/s,score[i][3]))
					else
						if pymolMode
							hits=0
						end
						push!(roc,(id1,id2,hits/s,score[i][3]))
					end
				end
			end
			return roc
		end
	end

	######################################################################
	# FUNCTION:		 getHmmLength
	######################################################################
	function getHmmLength(hmmFile::AbstractString)
		!isfile(hmmFile) && error("File not readable")
		for l in eachline(open(hmmFile))
			if l[1:4] == "LENG" || l[1:4]=="CLEN"
				LENG=parse(Int64,match(r"\d+",l).match)
				return LENG
			end
		end
		error("unable to get hmm length in $hmmFile")
	end

	######################################################################
	# FUNCTION:		 makeMarriedContactMap
	######################################################################
	## IDIOCY: THIS DOES THE SAME THING AS THE FUNCTION "interAlignDist"
	function makeMarriedContactMap(chain1::Chain,chain2::Chain;output::AbstractString="default")
		chain1map=chain1.mappedTo
		chain2map=chain2.mappedTo
		chain1map=="" && error("chain $(chain1.identifier) has no mapping")
		chain2map=="" && error("chain $(chain2.identifier) has no mapping")
		LENG1=getHmmLength(chain1map)
		LENG2=getHmmLength(chain2map)
		contactMap=-ones(Float64,LENG1+LENG2,LENG1+LENG2)
		# Protein 1
		for r1 in values(chain1.align)
			for r2 in values(chain1.align)
				r1.alignmentPos==-1 && error("incoherent alignment information")
				r2.alignmentPos==-1 && error("incoherent alignment information")
				contactMap[r1.alignmentPos,r2.alignmentPos]=residueDist(r1,r2)
			end
		end
		# Protein 2
		for r1 in values(chain2.align)
			for r2 in values(chain2.align)
				r1.alignmentPos==-1 && error("incoherent alignment information")
				r2.alignmentPos==-1 && error("incoherent alignment information")
				contactMap[r1.alignmentPos+LENG1,r2.alignmentPos+LENG1]=residueDist(r1,r2)
			end
		end
		# Protein 1 - Protein 2
		for r1 in values(chain1.align)
			for r2 in values(chain2.align)
				r1.alignmentPos==-1 && error("incoherent alignment information")
				r2.alignmentPos==-1 && error("incoherent alignment information")
				d=residueDist(r1,r2)
				contactMap[r1.alignmentPos,r2.alignmentPos+LENG1]=d
				contactMap[r2.alignmentPos+LENG1,r1.alignmentPos]=d
			end
		end
		if output!="default"
			fid=open(output,"w")
			for i=1:LENG1+LENG2
				for j=(i+1):LENG1+LENG2
					contactMap[i,j]<0.0 && continue
					@printf(fid,"%d\t%d\t%f\n",i,j,contactMap[i,j])
				end
			end
			close(fid)
		end
		return contactMap

	end

	######################################################################
	# FUNCTION:		 countContacts
	######################################################################
	function countContacts(chain::Chain;min_separation=5,cutoff=8.0)
		nums=sort([n for n in keys(chain.align)])
		contacts=0
		for i=1:length(nums)
			for j=(i+1):length(nums)
				id1=nums[i]
				id2=nums[j]
				if abs(id1-id2)<min_separation
					continue
				end
				res1=chain.align[id1]
				res2=chain.align[id2]
				if res1.alignmentPos<1 || res2.alignmentPos<1
					error("Something went horribly wrong here")
				end
				d=residueDist(res1,res2)
				if d<cutoff
					contacts=contacts+1
				end
			end
		end
		return contacts
	end

	######################################################################
	# FUNCTION:		 st2fa
	######################################################################
	function st2fa(iFile;oFile="default")
		!isfile(iFile) && error("File not found")
		if oFile=="default"
			oFile="$iFile.st2fa"
		end
		iFid=open(iFile,"r")
		oFid=open(oFile,"w")
		seqDict=Dict{AbstractString,AbstractString}()
		for line in eachline(iFid)
			line[1]=='#' && continue
			s=split(line)
			length(s)!=2 && continue
			if haskey(seqDict,s[1])
				seqDict[s[1]]*=s[2]
			else
				seqDict[s[1]]=s[2]
			end
		end
		close(iFid)
		for k in keys(seqDict)
			@printf(oFid,">%s \n%s\n",k,seqDict[k])
		end
		close(oFid)
		return oFile
	end

	######################################################################
	# FUNCTION:		 fasta2seqDict
	######################################################################
	function fasta2seqDict(iFile::AbstractString)
		!isfile(iFile) && error("File not found")
		iFid=open(iFile,"r")
		seqDict=Dict{AbstractString,AbstractString}()
		for head in eachline(iFid)
			( head[1]!='>' || eof(iFid) ) && error("Fasta not readable")
			seq=readline(iFid)
			haskey(seqDict,head) && error("Fasta contains several entries for $head")
			seqDict[head]=seq
		end
		close(iFid)
		return seqDict
	end

	######################################################################
	# FUNCTION:		 interactionSurface(chain1,chain2)
	######################################################################
	function interactionSurface(chain1,chain2;numbersOnly::Bool=true,alignedOnly::Bool=true,cutoff::Float64=8.0)
		if chain1.mappedTo=="" || chain2.mappedTo==""
			return interactionSurface_unmapped(chain1,chain2;numbersOnly=numbersOnly,cutoff=cutoff)
		end
		pdbPairs=0;
		alignmentPairs=0;
		if !numbersOnly
			iS=Array((AbstractString,AbstractString),0)
		end
		for r1 in values(chain1.residue)
			for r2 in values(chain2.residue)
				if residueDist(r1,r2)<cutoff
					pdbPairs+=1;
					if !numbersOnly && !alignedOnly
						push!(iS,(r1.identifier,r2.identifier))
					end
					if r1.alignmentPos>0 && r2.alignmentPos>0
						alignmentPairs+=1;
						if !numbersOnly && alignedOnly
							push!(iS,(r1.identifier,r2.identifier))
						end
					end

				end
			end
		end
		if numbersOnly
			return pdbPairs,alignmentPairs
		else
			return iS
		end

	end
	function interactionSurface_unmapped(chain1,chain2;numbersOnly::Bool=true,cutoff::Float64=8.0)
		pdbPairs=0;
		if !numbersOnly
			iS=Array((String,String),0)
		end
		for r1 in values(chain1.residue)
			for r2 in values(chain2.residue)
				if residueDist(r1,r2)<cutoff
					pdbPairs+=1;
					if !numbersOnly
						push!(iS,(r1.identifier,r2.identifier))
					end
				end
			end
		end
		if numbersOnly
			return pdbPairs
		else
			return iS
		end

	end

#</module>
end


text:
## Polynomial iterator
start(poly::Polynomial) = start(poly.poly)
next(poly::Polynomial, state) = next(poly.poly, state)
done(poly::Polynomial, state) = done(poly.poly, state)
length(poly::Polynomial) = length(poly.poly)
haskey(poly::Polynomial, key) = haskey(poly.poly, key)
keys(poly::Polynomial) = keys(poly.poly)
values(poly::Polynomial) = values(poly.poly)
getindex(poly::Polynomial, expo::Exponent) = poly.poly[expo]

## Exponent iterator
start(expo::Exponent) = start(expo.expo)
next(expo::Exponent, state) = next(expo.expo, state)
done(expo::Exponent, state) = done(expo.expo, state)
length(expo::Exponent) = length(expo.expo)
haskey(expo::Exponent, key) = haskey(expo.expo, key)
keys(expo::Exponent) = keys(expo.expo)
values(expo::Exponent) = values(expo.expo)
getindex(expo::Exponent, var::Variable) = expo.expo[var]

##################################
## Addition
##################################
function add!(p::Polynomial, p1::T) where T<:AbstractPolynomial
    return add!(p, convert(Polynomial, p1))
end

function add!(p::Polynomial, p1::Polynomial)
    for (cur_expo, λ) in p1
        λ != 0 || continue
        add_to_dict!(p.poly, cur_expo, λ)
    end
    p.degree.explvar = max(p.degree.explvar, p1.degree.explvar)
    p.degree.conjvar = max(p.degree.conjvar, p1.degree.conjvar)
end

function add(p1::Polynomial, p2::Polynomial)
    p = copy(p1)
    add!(p, p2)
    return p
end

function +(p1::T, p2::U) where T<:AbstractPolynomial where U<:AbstractPolynomial
    return add(convert(Polynomial, p1), convert(Polynomial, p2))
end
function +(p1::Number, p2::T) where T<:AbstractPolynomial
    return add(convert(Polynomial, p1), convert(Polynomial, p2))
end
function +(p1::T, p2::Number) where T<:AbstractPolynomial
    return add(convert(Polynomial, p1), convert(Polynomial, p2))
end

##################################
## Product
##################################
function product(p1::Polynomial, p2::Polynomial)
    p = Polynomial()
    for (expo1, λ1) in p1
        λ1 != 0 || continue
        for (expo2, λ2) in p2
            λ2 != 0 || continue
            expoprod = product(expo1, expo2)
            add_to_dict!(p.poly, expoprod, λ1 * λ2)
        end
    end
    p.degree.explvar = p1.degree.explvar+p2.degree.explvar
    p.degree.conjvar = p1.degree.conjvar+p2.degree.conjvar
    return p
end

function product(exp1::Exponent, exp2::Exponent)
    function add_expod!(expod, expod1)
        for (var, deg) in expod1
            if !haskey(expod, var)
                expod[var] = Degree(0,0)
            end
            expod[var].explvar += deg.explvar
            expod[var].conjvar += deg.conjvar
            if (expod[var].explvar, expod[var].conjvar) == (0,0)
                delete!(expod, var)
            end
        end
        return expod
    end

    expod = Dict{Variable, Degree}()
    add_expod!(expod, exp1.expo)
    add_expod!(expod, exp2.expo)
    return Exponent(expod)
end

function *(p1::T, p2::U) where T<:AbstractPolynomial where U<:AbstractPolynomial
    return product(convert(Polynomial, p1), convert(Polynomial, p2))
end
function *(p1::Number, p2::T) where T<:AbstractPolynomial
    return product(convert(Polynomial, p1), convert(Polynomial, p2))
end
function *(p1::T, p2::Number) where T<:AbstractPolynomial
    return product(convert(Polynomial, p1), convert(Polynomial, p2))
end

### Soustraction
function -(p::T) where T<:AbstractPolynomial
    return -1 * convert(Polynomial, p)
end

function -(p1::T, p2::U) where T<:AbstractPolynomial where U<:AbstractPolynomial
    return p1 + (-p2)
end
function -(p1::Number, p2::T) where T<:AbstractPolynomial
    return p1 + (-p2)
end
function -(p1::T, p2::Number) where T<:AbstractPolynomial
    return p1 + (-p2)
end

## Division
function divide(p1::Polynomial, p2::Polynomial)
    if length(p2) != 1
        error("/(::Polynomial, ::Polynomial): Only allowed for monomial divisor ($(length(p2))-monomial polynomial here).")
    end
    expo, λ = collect(p2)[1]
    if expo.degree != Degree(0,0)
        error("/(::Polynomial, ::Polynomial): Only allowed for constant divisor ($(expo.degree)-degree monomial here).")
    end
    if λ == 0
        error("/(::Polynomial, ::Polynomial): Only allowed for non null constant divisor.")
    end
    return p1 * (1/λ)
end

function /(p1::T, p2::U) where T<:AbstractPolynomial where U<:AbstractPolynomial
    return divide(convert(Polynomial, p1), convert(Polynomial, p2))
end
function /(p1::Number, p2::T) where T<:AbstractPolynomial
    return divide(convert(Polynomial, p1), convert(Polynomial, p2))
end
function /(p1::T, p2::Number) where T<:AbstractPolynomial
    return divide(convert(Polynomial, p1), convert(Polynomial, p2))
end

## Exponent
function powpoly(p::Polynomial, d::Int)
    if length(p) != 1
        error("^(::Polynomial, ::Int): Only supported for monomials ($(length(p))-monomial polynomial here).")
    end
    expo, λ = collect(p)[1]

    expod = Dict{Variable, Degree}()
    for (var, deg) in expo.expo
        expod[var] = Degree(deg.explvar*d, deg.conjvar*d)
    end
    return Polynomial(Dict{Exponent, Number}(Exponent(expod)=>λ^d))
end

function ^(p::T, d::Int) where T<:AbstractPolynomial
    return powpoly(convert(Polynomial, p), d)
end

ext: jl

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
import Pkg

try
    Pkg.add("Coverage")
catch e1
    @warn(
        string("Ignoring exception e1:"),
        e1,
        )
end

import Coverage

import RepoSnapshots

cd(RepoSnapshots.package_directory())

Coverage.Codecov.submit(Coverage.Codecov.process_folder())


text:
using Distributions
include("utils.jl")

"""
Given a 2D array of strings *sequences*, greedy
search for a motif of length *l*.
Return an array of indices corresponding to the
alignment found with the highest score.
"""
function motifsearch_greedy(sequences, l)
    # TODO: Check that l <= sequence length
    bestmotif = ones(Int32, 2)

    for s₁ = 1:length(sequences[1])-l+1
        for s₂ = 1:length(sequences[2])-l+1
            # find best partial alignment of first 2 sequences
            if score(sequences, [s₁, s₂], l) > score(sequences, bestmotif, l)
                bestmotif = [s₁, s₂]
            end
        end
    end
    s = copy(bestmotif)

    #iterate through the rest of the sequences
    for i = 3:length(sequences)
        # extend the default best motif
        push!(bestmotif, 1)
        for sᵢ = 1:length(sequences[i])-l+1
            if score(sequences, [s..., sᵢ], l) > score(sequences, bestmotif, l)
                # update the best motif
                bestmotif[i] = sᵢ
            end
        end
        push!(s, bestmotif[i])
    end
    return bestmotif
end

"""
Given a 2D array of strings *sequences*, randomized greedy
search for a motif of length *l*.
Return an array of indices corresponding to the
alignment found with the highest score.
"""
function motifsearch_random(sequences, l)
    # TODO: Check that l <= sequence length

    n = length(sequences[1])

    # generate an initial random alignment
    s = rand(1:n-l+1, length(sequences))

    # form profile P from s
    P = getprofile(align(sequences, s, l))

    bestscore = 0
    newscore = score(P)
    while newscore > bestscore
        bestscore = newscore
        for i = 1:length(sequences)
            # find a P-most probable l-mer from the ith sequence
            a = map(x -> probscore(x, P), [sequences[i][j:j+l-1] for j=1:n-l+1])
            # update the starting index for sᵢ
            s[i] = indmax(a)
        end
        # recalculate the profile
        P = getprofile(align(sequences, s, l))
        newscore = score(P)
    end

    return s
end

"""
Given a 2D array of strings *sequences*, randomized greedy
search for a motif of length *l*. *leniency* indicates how many steps
are taken after score stops increasing.
Return an array of indices corresponding to the alignment found with the
highest score.
"""
function motifsearch_gibbs(sequences, l, leniency)
    # TODO: Check that l <= sequence length

    n = length(sequences[1])

    # generate an initial random alignment
    s⃗ = rand(1:n-l+1, length(sequences))
    bests⃗ = copy(s⃗)

    bestscore = score(sequences, s⃗, l)
    currscore = bestscore
    # watch for convergence
    steps_since_inc = 0

    # do until convergence:
    while steps_since_inc < leniency
        # Randomly select one sequence, seqᵢ
        seqᵢ = rand(1:length(sequences))
        # Get the indices of the sequences left over
        inds = filter(x -> x != seqᵢ, 1:length(sequences))
        # Create profile from sequences not selected
        alignment = align(sequences[inds], s⃗[inds], l)
        profile = getprofile(alignment)
        # Calculate P-prob for each motif in seqᵢ
        a = map(x -> probscore(x, profile), [sequences[seqᵢ][j:j+l-1] for j=1:n-l+1])
        # if all zeros, set to 1 for uniform dist
        if all(a .== 0)
            a = ones(length(a))
        end
        p⃗ = a ./ sum(a)
        # Choose new starting position, sampled from probs found for seqᵢ
        c = Categorical(p⃗)
        s⃗[seqᵢ] = rand(c)

        # calculate new alignment score and check for convergence
        currscore = score(sequences, s⃗, l)
        if currscore <= bestscore
            # increment convergence check
            steps_since_inc += 1
        else
            # update bests⃗ if higher score found
            bestscore = currscore
            bests⃗ = copy(s⃗)
            steps_since_inc = 0
        end

    end

    return bests⃗
end

max_stars_repo_name: dillondaudert/JuliaMotifs

max_issues_repo_name: dillondaudert/JuliaMotifs

max_forks_repo_name: dillondaudert/JuliaMotifs


text:
using Zygote, OptimKit
using ..ITensorAutoHOOT
using ..ITensorNetworks
using ..ITensorAutoHOOT: batch_tensor_contraction
using ..ITensorNetworks:
  PEPS, inner_network, inner_networks, flatten, rayleigh_quotient, tree, Models

function loss_grad_wrap(peps::PEPS, Hs::Array, tensortype; kwargs...)
  function loss(peps::PEPS)
    peps_bra = addtags(linkinds, peps, "bra")
    peps_ket = addtags(linkinds, peps, "ket")
    sites = commoninds(peps_bra, peps_ket)
    peps_ket_ham = prime(sites, peps_ket)
    network_H = inner_networks(peps_bra, peps_ket, peps_ket_ham, Hs)
    network_inner = inner_network(peps_bra, peps_ket)
    network_list = vcat(network_H, [network_inner])
    variables = flatten([peps_bra, peps_ket, peps_ket_ham])
    inners = batch_tensor_contraction(tensortype, network_list, variables...; kwargs...)
    return rayleigh_quotient(inners)
  end
  loss_w_grad(peps::PEPS) = loss(peps), gradient(loss, peps)[1]
  return loss_w_grad
end

function loss_grad_wrap(peps::PEPS, Hs::Array{Models.LineMPO}, tensortype, tree; kwargs...)
  Hs_row = [H for H in Hs if H.coord[2] isa Colon]
  Hs_column = [H for H in Hs if H.coord[1] isa Colon]
  function loss(peps::PEPS)
    peps_bra = addtags(linkinds, peps, "bra")
    peps_ket = addtags(linkinds, peps, "ket")
    sites = commoninds(peps_bra, peps_ket)
    peps_ket_ham = prime(sites, peps_ket)
    # generate network
    network_list_row = inner_networks(peps_bra, peps_ket, peps_ket_ham, Hs_row, tree)
    network_list_column = inner_networks(peps_bra, peps_ket, peps_ket_ham, Hs_column, tree)
    network_inner = inner_network(peps_bra, peps_ket, tree)
    network_list = vcat(network_list_row, network_list_column, [network_inner])
    variables = flatten([peps_bra, peps_ket, peps_ket_ham])
    inners = batch_tensor_contraction(
      tensortype, network_list, variables...; optimize=false, kwargs...
    )
    return rayleigh_quotient(inners)
  end
  loss_w_grad(peps::PEPS) = loss(peps), gradient(loss, peps)[1]
  return loss_w_grad
end

function gradient_descent(
  peps::PEPS, Hs::Array, tensortype; stepsize::Float64, num_sweeps::Int, kwargs...
)
  loss_w_grad = loss_grad_wrap(peps, Hs, tensortype; kwargs...)
  return gradient_descent(peps, loss_w_grad; stepsize=stepsize, num_sweeps=num_sweeps)
end


text:
const AbstractTensor{N,T} = AbstractArray{T,N}
const TupleN{T,N} = NTuple{N,T}

struct Hankel{T,N,C,K,A<:AbstractArray} <: AbstractArray{T,N}
    image::A
    channel_size::NTuple{C,Int}
    kernel_size::NTuple{K,Int}
    @doc raw"""
        Hankel(image, kernel_size, channel_size)

    Constructs a Hankel array `A` from the underlying `image` data, such that:

    ```julia
    A[c,j,k,b] = image[c,j+k-1,b]
    ```

    where `c,j,k,b` are multi-indices for the channel, kernel, image, and batch
    dimensions. So `j` is the kernel entry, and `k` the kernel position in the image.
    """
    function Hankel(image::AbstractArray, channel_size::NTuple{C,Int}, kernel_size::NTuple{K,Int}) where {C,K}
        input_size = size(image)[(C + 1):(C + K)]
        batch_size = size(image)[(C + K + 1):end]
        @assert size(image) == (channel_size..., input_size..., batch_size...)
        @assert all(input_size .≥ kernel_size)
        return new{eltype(image), ndims(image) + K, C, K, typeof(image)}(image, channel_size, kernel_size)
    end
end

function Hankel(image::AbstractArray, C::Int, kernel_size::NTuple{K,Int}) where {K}
    channel_size = size(v)[1:C]
    return Hankel(image, channel_size, kernel_size)
end

Base.IndexStyle(::Type{<:Hankel}) = IndexCartesian()
Base.size(A::Hankel) = (channel_size(A)..., kernel_size(A)..., output_size(A)..., batch_size(A)...)

image_size(A::Hankel) = size(A.image)
channel_size(A::Hankel) = A.channel_size
kernel_size(A::Hankel) = A.kernel_size
input_size(A::Hankel) = image_size(A)[(channel_ndims(A) + 1):(channel_ndims(A) + kernel_ndims(A))]
output_size(A::Hankel) = input_size(A) .- kernel_size(A) .+ 1
batch_size(A::Hankel) = image_size(A)[(channel_ndims(A) + kernel_ndims(A) + 1):end]

channel_ndims(A::Hankel) = length(channel_size(A))
kernel_ndims(A::Hankel) = length(kernel_size(A))

@inline function Base.getindex(A::Hankel{T,N}, I::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, I...)
    @inbounds begin
        c = I[1:channel_ndims(A)]
        j = I[(channel_ndims(A) + 1):(channel_ndims(A) + kernel_ndims(A))]
        k = I[(channel_ndims(A) + kernel_ndims(A) + 1):(channel_ndims(A) + 2kernel_ndims(A))]
        b = I[(channel_ndims(A) + 2kernel_ndims(A) + 1):end]
        i = j .+ k .- 1
        return A.image[CartesianIndex(c), CartesianIndex(i), CartesianIndex(b)]
    end
end

@doc raw"""
    hankel(v, cdims, ksz)

Creates a Hankel array from `v` with channel dimensions `cdims` and kernel size `ksz`.
The returned array `A` satisfies:

```math
A_{i,j,k,n} = v_{i,j+k-1,n}
```

where `i,j,k,n` are multi-indices, with `i` traversing channels
and `j` traversing the kernel.
```
"""
function hankel(img::AbstractArray, channel_size::NTuple{C,Int}, kernel_size::NTuple{K,Int}) where {C,K}
    input_size = size(img)[(C + 1):(C + K)]
    batch_size = size(img)[(C + K + 1):end]
    @assert size(img) == (channel_size..., input_size..., batch_size...)
    output_size = input_size .- kernel_size .+ 1
    A = similar(img, eltype(img), channel_size..., kernel_size..., output_size..., batch_size...)
    vflat = reshape(img, prod(channel_size), input_size..., prod(batch_size))
    Aflat = reshape(A, prod(channel_size), kernel_size..., output_size..., prod(batch_size))
    flat_hankel!(Aflat, vflat)
    return A
end

function hankel(v::AbstractArray, C::Int, kernel_size::TupleN{Int})
    channel_size = size(v)[1:C]
    return hankel(v, channel_size, kernel_size)
end

# for 1 channel dimension and 1 batch dimension
flat_hankel!(A::AbstractTensor{4}, v::AbstractTensor{3}) =
    @tullio A[c,j1,k1,b] = v[c, j1 + k1 - 1, b]
flat_hankel!(A::AbstractTensor{6}, v::AbstractTensor{4}) =
    @tullio A[c,j1,j2,k1,k2,b] = v[c, j1 + k1 - 1, j2 + k2 - 1, b]
flat_hankel!(A::AbstractTensor{8}, v::AbstractTensor{5}) =
    @tullio A[c,j1,j2,j3,k1,k2,k3,b] = v[c, j1 + k1 - 1, j2 + k2 - 1, j3 + k3 - 1, b]


text:
#=validate.jl - replay the game from start to the current move, at every step making
sure that the rules were kept. Print 0 is the game was played cleanly so far and the move ID
of the violating move if it wasn’t. Accepts 1 command line argument,<filename> => database
=#

#=
Only captured targets can be dropped
Dropped cannot capture a piece
It is not promoted instantly
Pawn, knight, lance cannot be dropped on the furthest rank
Pawn cannot be dropped on the same column (y) as another unpromoted pawn
Pawn cannot checkmate only checks
=#
module validateMod

include("board.jl")
include("dParse.jl")
gameType = "standard" #TO DO: don't hard code this
  using BM
  using SQLite
  board = (BM.startGame(gameType)).state
  database = ARGS[1] #/path/to/database/file {string}
  db = SQLite.DB(database) #Opens the database gamefile

  #returns True if move is Valid, False otherwise
  function moveValidate(unit, moveType, gameType, team, sourcex, sourcey, targetx, targety, targetx2, targety2)

    #edge 0, sourcex or source y are null but it is not drop
    if ( (sourcex == -1 || sourcey == -1) && moveType != "drop")
      return false
    end

    #edge 1, source and target are the same
    if (sourcex == targetx) && (sourcey == targety)
      return false
    end

    #edge 2, source or target are out of bounds
    if (gameType == "standard") #normalShogi
      if  (~(sourcex >= 1 && sourcex <= 9) || ~(sourcey >= 1 && sourcey <= 9) || ~(targetx >= 1 && targetx <= 9) || ~(targety >= 1 && targety <= 9))
        return false
      end
    elseif (gameType == "minishogi") #miniShogi
      if  (~(sourcex >= 1 && sourcex <= 5) || ~(sourcey >= 1 && sourcey <= 5) || ~(targetx >= 1 && targetx <= 5) || ~(targety >= 1 && targety <= 5))
        return false
      end
    elseif (gameType == "chu") #chuShogi
      if  (~(sourcex >= 1 && sourcex <= 12) || ~(sourcey >= 1 && sourcey <= 12) || ~(targetx >= 1 && targetx <= 12) || ~(targety >= 1 && targety <= 12))
        return false
      end
    end

    #drop case
    if (moveType == "drop")

      #Chu Shogi has no drops
      if (moveType == "chu")
        return false
      end

      #drop case 1: Dropped cannot capture
      if (targetx,targety) in keys(board)
        return false
      end

      #drop case 2: Pawn, Knight, Lance cannot be dropped on the furthest rank
      if (unit == "pawn" || unit == "knight" || unit == "lance")
        if (team == 'b') #team black
          if (targetx == 1)
            return false
          end
        elseif (team =='w') #team white
          if (targetx == 9)
            return false
          end
        end
      end

      #drop case3: It is not promoted
      if (unit == "promotedbishop" || unit == "promotedlance" || unit == "promotedknight" || unit == "promotedpawn" || unit == "promotedrook" || unit == "promotedsilver general")
        return false
      end

      #drop case4: Pawns cannot be dropped on the same column as another unpromoted pawn
      if (unit == "pawn")
        for x in 1:9
          if (board[(x,targety)][1] == "pawn") #there is a pawn
            return false
          end
        end
      end

      #drop case5: Pawns cannot checkmate #missing

      #drop case6: Capture targets can be dropped #missing

    end

    #case 1 bishop
    if unit == "bishop"
      return bishopValidate(team,sourcex,sourcey,targetx,targety)

    #case 1.2 promoted bishop or dragon Horse
    elseif unit == "promotedbishop"
        return pBishopValidate(team,sourcex,sourcey,targetx,targety)
      end

    #case 2 gold general
    elseif unit == "goldgeneral"
      return goldGeneralValidate(team,sourcex,sourcey,targetx,targety)

    #case 3 king
    elseif unit == "king"
      return kingValidate(team,sourcex,sourcey,targetx,targety)

    #case 4 lance
    elseif unit == "lance"
      return lanceValidate(team,sourcex,sourcey,targetx,targety)

    #case 4.2 promoted lance
    elseif unit == "promotedlance"
      #promotedLance moves same as gold general
      return goldGeneralValidate(team,sourcex,sourcey,targetx,targety)

    #case 5 knight or Kirin
    elseif unit == "knight"
      return knightValidate(team,sourcex,sourcey,targetx,targety)

    elseif unit == "kirin" && gameType == "chu"
      return kirinValidate(team,sourcex,sourcey,targetx,targety)

    #case 5.2 promoted knight
    elseif unit == "promotedknight"
      #promotedKnight moves same as gold general
      return goldGeneralValidate(team,sourcex,sourcey,targetx,targety)

    #case 6 pawn
    elseif unit == "pawn"
      return pawnValidate(team,sourcex,sourcey,targetx,targety)

    #case 6.2 promoted pawn
    elseif unit == "promotedpawn"
      #promotedPawn moves same as gold general
      return goldGeneralValidate(team,sourcex,sourcey,targetx,targety)

    #case 7 rook
    elseif unit == "rook"
      return rookValidate(team,sourcex,sourcey,targetx,targety)

    #case 7.2 promoted rook or dragon king
    elseif unit == "promotedrook"
      return pRookValidate(team,sourcex,sourcey,targetx,targety)

    #case 8 silver general
    elseif unit == "silvergeneral"
      return silverGeneralValidate(team,sourcex,sourcey,targetx,targety)

    #case 8.2 promoted silver general or veritcal mover
    elseif unit == "promotedsilver general"
      #promoted silver general moves same as gold general
      return goldGeneralValidate(team,sourcex,sourcey,targetx,targety)

    #CHU SHOGI
    #case 9 Reverse Chariot promotes to whale
    elseif unit == "reversechariot" && gameType == "chu"
      return reverseChariotValidate(team,sourcex,sourcey,targetx,targety)

    #case 9.2 whale
    elseif unit == "whale" && gameType == "chu"
      return whaleValidate(team,sourcex,sourcey,targetx,targety)

    #case 10 Copper general promotes to side mover
    elseif unit == "coppergeneral" && gameType == "chu"
      return copperGeneralValidate(team,sourcex,sourcey,targetx,targety)

    #case 11 Dragon King promotes to soaring eagle
    elseif unit == "dragonking" && gameType == "chu"
      return dragonKingValidate(team,sourcex,sourcey,targetx,targety)

    #case 11.2 Soaring eagle
    elseif unit == "soaringeagle" && gameType == "chu"
      return soaringEagleValidate(team,sourcex,sourcey,targetx,targety,targetx2,targety2)

    #case 12 Drunk Elephant promotes to prince
    elseif unit == "drunkelephant" && gameType == "chu"
      return drunkElephantValidate(team,sourcex,sourcey,targetx,targety)

    #case 12.2 prince
    elseif unit == "prince" && gameType == "chu"
      #move like king
      return kingValidate(team,sourcex,sourcey,targetx,targety)

    #case 13 ferocious leopard promotes bishop
    elseif unit == "ferociousleopard" && gameType == "chu"
      return ferociousLeopardValidate(team,sourcex,sourcey,targetx,targety)

    #case 14 Dragon Horse promotes to horned falcon
    elseif unit == "dragonhorse" && gameType == "chu"
      return dragonHorseValidate(team,sourcex,sourcey,targetx,targety)

    #case 14.2 horned falcon
    elseif unit == "hornedfalcon" && gameType == "chu"
      return hornedFalconValidate(team,sourcex,sourcey,targetx,targety,targetx2,targety2)

    #case 15 Lion
    elseif unit == "lion" && gameType == "chu"
      return lionValidate(team,sourcex,sourcey,targetx,targety,targetx2,targety2)

    #case 16 Side Mover promotes to free boar
    elseif unit == "sidemover" && gameType == "chu"
      return sideMoverValidate(team,sourcex,sourcey,targetx,targety)

    #case 16.2 free boar
    elseif unit == "freeboar" && gameType == "chu"
      return freeBoarValidate(team,sourcex,sourcey,targetx,targety)

    #case 17 Go-between promotes to drunk Elephant
    elseif unit == "gobetween" && gameType == "chu"
      return goBetweenValidate(team,sourcex,sourcey,targetx,targety)

    #case 18 Blind tiger promotes to flying stag
    elseif unit == "blindtiger" && gameType == "chu"
      return blindTigerValidate(team,sourcex,sourcey,targetx,targety)

    #case 18.2 flying stag
    elseif unit == "flyingstag" && gameType == "chu"
      return flyingStagValidate(team,sourcex,sourcey,targetx,targety)

    #case 19 Queen
    elseif unit == "queen" && gameType == "chu"
      return queenValidate(team,sourcex,sourcey,targetx,targety)

    #case 20 Vertical Mover promotes to flying ox
    elseif unit == "verticalmover" && gameType == "chu"
      return verticalMoverValidate(team,sourcex,sourcey,targetx,targety)

    #case 20.2
    elseif unit == "flyingox" && gameType == "chu"
      return flyingOxValidate(team,sourcex,sourcey,targetx,targety)

    #case 21 Phoenix promotes to queen
    elseif unit == "phoenix" && gameType == "chu"
      return phoenixValidate(team,sourcex,sourcey,targetx,targety)

    else
      return false #no valid unit
    end
  end

  #moveUpValidate <=== Return true if unit can move veritcally upwards to target ===>
  function moveUpValidate(sourcex, sourcey, targetx, targety)
    if (sourcey != targety)
      return false
    end

    unitCheck = abs(sourcex - targetx) - 1
    x = sourcex - 1
    for unit in 1:unitCheck
      if !((x,targety) in keys(board)) #Empty
        x = x - 1
      else
        return false
      end
    end
    return true
  end #moveUpValidate

  #moveDownValidate <=== Return true if unit can move veritcally downwards to target ===>
  function moveDownValidate(sourcex, sourcey, targetx, targety)
    if (sourcey != targety)
      return false
    end

    unitCheck = abs(sourcex - targetx) - 1
    x = sourcex + 1
    for unit in 1:unitCheck
      if !((x,targety) in keys(board)) #isEmpty
        x = x + 1
      else
        return false
      end
    end

    return true
  end #moveDownValidate

  #moveLeftValidate <=== Return true if unit can move horizontally left to target ===>
  function moveLeftValidate(sourcex, sourcey, targetx, targety)
    if (sourcex != targetx)
      return false
    end

    unitCheck = abs(sourcey - targety) - 1
    y = sourcey - 1
    for unit in 1:unitCheck
      if !((targetx,y) in keys(board)) #isEmpty
        y = y - 1
      else
        return false
      end
    end

    return true
  end #moveLeftValidate

  #moveRightValidate <=== Return true if unit can move horizontally Right to target ===>
  function moveRightValidate(sourcex, sourcey, targetx, targety)
    if (sourcex != targetx)
      return false
    end

    unitCheck = abs(sourcey - targety) - 1
    y = sourcey + 1
    for unit in 1:unitCheck
      if !((targetx,y) in keys(board))
        y = y + 1
      else
        return false
      end
    end

    return true
  end #moveLeftValidate

  #moveOrthogonalValidate <=== Return true if unit can NESW directions ===>
  function moveOrthogonalValidate(sourcex,sourcey,targetx,targety)
    #checks for moves like rook
    if (sourcex == targetx) #horizontal
      if (sourcey < targety) #move right
        if (moveRightValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      else #move left
        if (moveLeftValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      end
    elseif (sourcey == targety) #vertical
      if (sourcex < targetx) #move down
        if (moveDownValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      else #move up
        if (moveUpValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      end
    end

    return false
  end

  #moveAdjacentValidate <=== Return true if unit can move 1 UNIT adjacent ===>
  function moveAdjacentValidate(sourcex,sourcey,targetx,targety)
    if ((abs(sourcex - targetx) == 1) || (abs(sourcex - targetx) == 0))
      if ((abs(sourcey - targety) == 1) || (abs(sourcey - targety) == 0))
        return true
      end
    end

    return false
  end

  #diagonalUpLeft <=== Return true if unit can move diagonally upleft to target ===>
  function diagonalUpLeftValidate(sourcex, sourcey, targetx, targety)
    if (sourcex - targetx) == (sourcey - targety)
      if (sourcex > targetx) #moving diagonal up left
        unitCheck == abs(sourcex - targetx) - 1
        x = sourcex - 1
        y = sourcey - 1
        for unit in 1:unitCheck
          if !((x,y) in keys(board)) #isEmpty
            x = x - 1
            y = y - 1
          else
            return false
          end
        end
        return true
      end
    end

    return false
  end #diagonalUpLeftValidate

  #digonalUpRight <=== Return true if unit can move diagonally upleft to target ===>
  function diagonalUpRightValidate(sourcex, sourcey, targetx, targety)
    if (sourcex - targetx) == (targety - sourcey)
      if (sourcex > targetx) #moving diagonal up right
        unitCheck == abs(sourcex - targetx) - 1
        x = sourcex - 1
        y = sourcey + 1
        for unit in 1:unitCheck
          if !((x,y) in keys(board)) #isEmpty
            x = x - 1
            y = y + 1
          else
            return false
          end
        end
        return true
      end
    end

    return false
  end #diagonalUpRightValidate

  #diagonalDownLeft <=== Return true if unit can move diagonally downleft to target ===>
  function diagonalDownLeftValidate(sourcex, sourcey, targetx, targety)
    if (sourcex - targetx) == (targety - sourcey)
      if (sourcex < targetx) #moving diagonal down left
        unitCheck == abs(sourcex - targetx) - 1
        x = sourcex + 1
        y = sourcey - 1
        for unit in 1:unitCheck
          if !((x,y) in keys(board)) #isEmpty
            x = x + 1
            y = y - 1
          else
            return false
          end
        end
        return true
      end
    end

    return false
  end #diagonalDownLeftValidate

  #diagonalDownRight <=== Return true if unit can move diagonally downright to target ===>
  function diagonalDownRightValidate(sourcex, sourcey, targetx, targety)
    if (sourcex - targetx) == (sourcey - targety)
      if (sourcex < targetx) #moving diagonal down left
        unitCheck == abs(sourcex - targetx) - 1
        x = sourcex + 1
        y = sourcey + 1
        for unit in 1:unitCheck
          if !((x, y) in keys(board)) #isEmpty
            x = x + 1
            y = y + 1
          else
            return false
          end
        end
        return true
      end
    end

    return false
  end #diagonalDownRightValidate

  #moveDiagonalValidate <=== Return true if unit can move diagonally===>
  function moveDiagonalValidate(sourcex, sourcey, targetx, targety)
    #moves like bishop
    if (abs(sourcex - targetx) == abs(sourcey - targety)) #moves diagonally
      if (sourcex < targetx) #downwards
        if (sourcey < targety) #diagonalDownRight
          if (diagonalDownRightValidate(sourcex,sourcey,targetx,targety) == true)
            return true
          end
        else #diagonalDownLeft
          if (diagonalDownLeftValidate(sourcex,sourcey,targetx,targety) == true)
            return true
          end
        end
      else #upwards
        if (sourcey < targety) #diagonalUpRightValidate
          if (diagonalUpRightValidate(sourcex,sourcey,targetx,targety) == true)
            return true
          end
        else #diagonalUpLeft
          if (diagonalUpLeftValidate(sourcex,sourcey,targetx,targety) == true)
            return true
          end
        end
      end
    end

    return false
  end

  #case 1 #moves diagonally by any
  function bishopValidate(team, sourcex, sourcey, targetx, targety)
    if (moveDiagonalValidate(sourcex,sourcey,targetx,targety) == true)
      return true
    end

    #Not moving diagonally
    return false
  end #function bishopValidate end

  #case 1.2 moves like king and bishop
  function pBishopValidate(team,sourcex,sourcey,targetx,targety)

    #checks if it moves like king
    if ((abs(sourcex - targetx) == 1) || (abs(sourcex - targetx) == 0))
      if ((abs(sourcey - targety) == 1) || (abs(sourcey - targety) == 0))
        return true
      end
    end

    #checks moves like normal bishop
    if (moveDiagonalValidate(sourcex,sourcey,targetx,targety) == true)
      return true
    end

    return false #does not move like king or bishop
  end #function pBishopValidate end

  # case 2 gold general #moves like king but no diagonal backwards
  function goldGeneralValidate(team,sourcex,sourcey,targetx,targety)

    #check if it moves like king
    if ((abs(sourcex - targetx) == 1) || (abs(sourcex - targetx) == 0))
      if (abs(sourcey - targety) == 1) || (abs(sourcey - targety) == 0)
        if (team == 'b') #team black
          if (sourcex - targetx == -1) && (sourcey - targety == 1) #bottom left
            return false
          elseif (sourcex - targetx == -1) && (sourcey - targety == -1) #bottom right
            return false
          end
        elseif (team == 'w') #team white
          if (sourcex - targetx == 1) && (sourcey - targety == 1)#top left
            return false
          elseif (sourcex - targetx == 1) && (sourcey - targety == -1) #top right
            return false
          end
        else #invalid team
          return false
        end
        #when it is in either team and not in one of the two units
        return true
      end
    end

    return false #does not move like goldGeneral
  end #goldGeneralValidate

  #case 3 king #moves any by 1
  function kingValidate(team,sourcex,sourcey,targetx,targety)
    if (moveAdjacentValidate(sourcex,sourcey,targetx,targety) == true)
      return true
    end

    return false
  end #kingValidate end

  #case 4 lance #moves forward by any
  #hopping check
  function lanceValidate(team,sourcex,sourcey,targetx,targety)
    if (team == 'b') #team black
      if (targetx < sourcex) && (targety == sourcey)
        if (moveUpValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      end
    elseif (team == 'w') #team white
      if (targetx > sourcex) && (targety == sourcey)
        if (moveDownValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      end
    else
      return false #invalid team
    end

    return false
  end #lanceValidate end

  #case 4.2 promoted lance #moves same as gold general

  #case 5 knight #moves like L
  function knightValidate(team,sourcex,sourcey,targetx,targety)
    if (team == 'b') #team black
      if (sourcex - targetx == 2)
        if (abs(targety - sourcey) == 1)
          return true
        end
      end
    elseif (team == 'w') #team white
      if (sourcex - targetx == -2)
        if (abs(targety - sourcey) == 1)
          return trues
        end
      end
    else #invalid team
      return false
    end

    return false
  end #knigthValidate

  #case 5.2 promoted knight #moves same as gold general

  #case 6 pawn
  function pawnValidate(team,sourcex,sourcey,targetx,targety)
    if (team == 'b') #team black
      if (targetx == sourcex - 1) && (targety == sourcey)
        return true
      end
    elseif (team == 'w') #team white
      if (targetx == sourcex + 1) && (targety == sourcey)
        return true
      end
    else
      return false #invalid team
    end

    return false
  end #pawnValidate end

  #case 6.2 promotedPawn #moves same as gold general

  #case 7 rook
  function rookValidate(team,sourcex,sourcey,targetx,targety)
    #moving horizontally
    if (moveOrthogonalValidate(sourcex,sourcey,targetx,targety) == true)
      return true
    end

    return false
  end #rookValidate end

  #case 7.2 promotedRook #moves like king or normal rook
  function pRookValidate(team,sourcex,sourcey,targetx,targety)
    #checks if it moves like king
    if ((abs(sourcex - targetx) == 1) || (abs(sourcex - targetx) == 0))
      if ((abs(sourcey - targety) == 1) || (abs(sourcey - targety) == 0))
        return true
      end
    end

    #checks if it moves like normal rook
    #moving horizontally
    if (sourcex == targetx) && (sourcey != targety)
      #horizontal left
      if (sourcey > targety)
        if (moveLeftValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      #horizontal right
      else #(sourcey < targety)
        if (moveRightValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      end

    #moving vertically
    elseif (sourcey == targety) && (sourcex != targetx)
      #vertical up
      if (sourcex > targetx)
        if (moveUpValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      #vertical down
      else #(sourcex < targetx)
        if (moveDownValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      end
    end

    return false
  end #pRookValidate end

  #case 8 silverGeneral
  function silverGeneralValidate(team,sourcex,sourcey,targetx,targety)
    res = false
    if (team == 'b') #team black
      if (abs(sourcex - targetx) == 1)
        if (abs(sourcey - targety) == 1 ) #every diagonal corner
          res = true
        elseif (sourcey == targety) && (sourcex == targetx + 1) #up 1
          res = true
        end
      end
    elseif (team == 'w') #team white
      #println("white")
      if (abs(sourcex - targetx) == 1)
        if (abs(sourcey - targety) == 1) #every diagonal corner
          res = true
        elseif (sourcey == targety) && (sourcex == targetx - 1) #down 1
          res = true
        end
      end
    end

    return res
  end #silverGeneralValidate end

  #case 8.2 pSilverGeneral moves same as gold general

  #case 9 reverseChariotValidate
  function reverseChariotValidate(team,sourcex,sourcey,targetx,targety)
    if (sourcey == targety)
      if (sourcex > targetx)
        if (moveUpValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      else #sourcex < targety
        if (moveDownValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      end
    end
    return false
  end #reverseChariotValidate end

  #case 9.2 whale
  function whaleValidate(team,sourcex,sourcey,targetx,targety)

    if (sourcey == targety)
      #moving up
      if (sourcex > targetx)
        if (moveUpValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      #moving down
      else
        if (moveDownValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      end
    elseif (abs(sourcex - targetx) == abs(sourcey - targety)) #diagonal
      if (team == 'b') && (sourcex < targetx)
        if (sourcey < targety) #diagonalDownRight
          if (diagonalDownRightValidate(sourcex,sourcey,targetx,targety) == true)
            return true
          end
        else #diagonalDownLeft
          if (diagonalDownLeftValidate(sourcex,sourcey,targetx,targety) == true)
            return true
          end
        end
      elseif (team == 'w') && (sourcex > targetx)
        if (sourcey < targety) #diagonalUpRight
          if (diagonalUpRightValidate(sourcex,sourcey,targetx,targety) == true)
            return true
          end
        else #diagonalUpLeft
          if (diagonalUpLeftValidate(sourcex,sourcey,targetx,targety) == true)
            return true
          end
        end
      end
    end

    return false
  end #whaleValidate end

  #case 10 copperGeneralValidate
  function copperGeneralValidate(team,sourcex,sourcey,targetx,targety)
    if (team == 'b')
      if (targetx == sourcex - 1)
        if ((sourcey == targety) || (abs(sourcey - targety) == 1)) #forward any direction by 1
          return true
        end
      elseif (targetx == sourcex + 1)
        if (sourcey == targety) #back 1
          return true
        end
      end
    elseif (team == 'w')
      if (targetx == sourcex + 1)
        if ((sourcey == targety) || (abs(sourcey - targety) == 1)) #forward any direction by 1
          return true
        end
      elseif (targetx == sourcex - 1)
        if (sourcey == targety)
          return true
        end
      end
    end

    return false
  end #copperGeneralValidate end

  #case 10.2 sideMoverValidate
  function sideMoverValidate(team,sourcex,sourcey,targetx,targety)
    if (sourcex == targetx)
      if (sourcey > targety) #moveLeftValidate
        if (moveLeftValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      else #moveRightValidate
        if (moveRightValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      end
    elseif (sourcey == targety)
      if (abs(sourcex - targetx) == 1)
        return true
      end
    end

    return false
  end #sideMoverValidate end

  #case 11 dragonKingValidate #moves like king and rook
  function dragonKingValidate(team,sourcex,sourcey,targetx,targety)
    #checks for moves like king
    if (moveAdjacentValidate(sourcex,sourcey,targetx,targety) == true)
      return true
    end

    #checks for moves like rook
    if (moveOrthogonalValidate(sourcex,sourcey,targetx,targety) == true)
      return true
    end

    return false
  end #dragonKingValidate end

  #case 11.2 Soaring eagle
  function soaringEagleValidate(team,sourcex,sourcey,targetx,targety,targetx2,targety2)
    #range move
    if (isnull(targetx2) && isnull(targety2))

      #moves like rook
      if (sourcex == targetx) && (sourcey != targety) #moving horizontally
        if (sourcey > targety) #horizontal left
          if (moveLeftValidate(sourcex,sourcey,targetx,targety) == true)
            return true
          end
        else #horizontal right
          if (moveRightValidate(sourcex,sourcey,targetx,targety) == true)
            return true
          end
        end
      elseif (sourcey == targety) && (sourcex != targetx) #moving vertically
        if (sourcex > targetx) #vertical up
          if (moveUpValidate(sourcex,sourcey,targetx,targety) == true)
            return true
          end
        else #vertical down
          if (moveDownValidate(sourcex,sourcey,targetx,targety) == true)
            return true
          end
        end
      end

      #diagonal down
      if (team == 'b')
        if (sourcex < targetx) #move down
          if (sourcey < targety) #move right
            if (diagonalDownRightValidate(sourcex,sourcey,targetx,targety) == true)
              return true
            end
          else #move left
            if (diagonalDownLeftValidate(sourcex,sourcey,targetx,targety) == true)
              return true
            end
          end
        else #move up diagonal 1 or 2
          if (sourcex - targetx == 1) && (abs(sourcey - targety) == 1)
            return true
          elseif (sourcex - targetx == 2) && (abs(sourcey - targety) == 2)
            return true
          end
        end

      elseif (team == 'w')
        if (sourcex > targetx) #move up
          if (sourcey < targety) #move right
            if (diagonalUpRightValidate(sourcex,sourcey,targetx,targety) == true)
              return true
            end
          else #move left
            if (diagonalUpLeftValidate(sourcex,sourcey,targetx,targety) == true)
              return true
            end
          end
        else #move up diagonal 1 or 2
          if (sourcex - targetx == -1) && (abs(sourcey - targety) == 1)
            return true
          elseif (sourcex - targetx == -2) && (abs(sourcey - targety) == 2)
            return true
          end
        end
      end

    else #lion move

      #case 1 capture without movin
      if (sourcex == targetx2 && sourcey == targety2)
        if (team == 'b')
          if ( (sourcex - targetx == 1) && (abs(sourcey - targety) == 1) )
            return true
          end
        elseif (team == 'w')
          if ( (targetx - sourcex == 1) && (abs(sourcey - targety) == 1) )
            return true
          end
        end
      end

      #case 2 eat 1 hop another
      if (team =='b')
        if (sourcex - targetx2 == 2) && (sourcex - targetx == 1)
          if (abs(sourcex - targetx) == 1) && (abs(sourcex - targetx2) == 2)
            return true
          end
        end
      elseif (team == 'w')
        if (targetx2 - sourcex == 2) && (targetx - sourcex == 1)
          if (abs(sourcex - targetx) == 1) && (abs(sourcex - targetx2) == 2)
            return true
          end
        end
      end
    end

    return false
  end #soaringEagleValidate end

  #case 12 drunkElephantValidate
  function drunkElephantValidate(team,sourcex,sourcey,targetx,targety)
    if ((abs(sourcex - targetx) == 1) || (abs(sourcex - targetx) == 0))
      if ((abs(sourcey - targety) == 1) || (abs(sourcey - targety) == 0))
        if (team == 'b')
          if (sourcey == targety) && (targetx == sourcex + 1) #moveback
            return false
          else
            return true
          end
        elseif (team == 'w')
          if (sourcey == targety) && (targetx == sourcex - 1) #moveUp
            return false
          else
            return true
          end
        end
      end
    end

    return false
  end #drunkElephantValidate

  #case 12.2 princeValidate moves like king

  #case 13 ferociousLeopardValidate
  function ferociousLeopardValidate(team,sourcex,sourcey,targetx,targety)
    if ((abs(sourcex - targetx) == 1) || (abs(sourcex - targetx) == 0))
      if ((abs(sourcey - targety) == 1) || (abs(sourcey - targety) == 0))
        if (sourcex != targetx)
          return true
        end
      end
    end

    return false
  end #ferociousLeopardValidate end

  #case 13.2 bishop

  #case 14 dragonHorseValidate
  function dragonHorseValidate(team,sourcex,sourcey,targetx,targety)
    #check for moves like king
    if (moveAdjacentValidate(sourcex,sourcey,targetx,targety) == true)
      return true
    end

    #check for moves like bishop
    if (moveDiagonalValidate(sourcex,sourcey, targetx, targety) == true)
      return true
    end

    return false
  end #dragonHorseValidate end

  #case 14.2 hornedFalconValidate
  function hornedFalconValidate(team,sourcex,sourcey,targetx,targety,targetx2,targety2)

    #range moves
    if isnull(targetx2) && isnull(targety2)

      #check if moves like bishop
      if (moveDiagonalValidate(sourcex,sourcey,targetx,targety) == true)
        return true
      end

      #Check for vertical down moves and jump 2 forward
      if (sourcey == targety)
        if (team == 'b')
          if (sourcex < targetx)
            if (moveDownValidate(sourcex,sourcey,targetx,targety) == true)
              return true
            end
          else #jump 2 up
            if (sourcex - targetx == 2)
              return true
            end
          end
        elseif (team == 'w')
          if (sourcex > targetx)
            if (moveUpValidate(sourcex,sourcey,targetx,targety) == true)
              return true
            end
          else #jump 2 down
            if (targetx - sourcex == 2)
              return true
            end
          end
        end
      end

    else #lion moves #where targetx2 and targety2 is NOT null

      #case 2 eats 1 forward without moving
      #case 3 eat 1 forward and then jump 1 forward # not your own team

      if (sourcey == targety)
        if (team == 'b')
          #case 2 eats 1 forward without moving
          if (sourcex == targetx2 && sourcey == targety2)
            if (sourcex == targetx + 1)
              return true
            end
          #case 3 eat 1 forward and then jump 1  forward
          elseif (sourcex == targetx2 + 2)
            if (sourcex == targetx + 1)
              #check if it has an enemy here
              if ((targetx,targety) in keys(board)) && (board[targetx,targety].team == 'w')
                return true
              end
            end
          end
        elseif (team == 'w')
          #case 2 eats 1 forward without moving
          if (sourcex == targetx2 && sourcey == targety2)
            if (sourcex == targetx - 1)
              return true
            end
          #case 3 eat 1 forward and then jump 1  forward
          elseif (sourcex == targetx2 - 2)
            if (sourcex == targetx - 1)
              #check if it has an enemy here
              if ((targetx,targety) in keys(board)) && (board[targetx,targety].team == 'b')
                return true
              end
            end
          end
        end
      end
    end

    return false
  end #hornedFalconValidate end

  #case 15 lionValidate
  function lionValidate(team,sourcex,sourcey,targetx,targety,targetx2,targety2)

    flag2 = false
    if (targetx2 != -1 || targety2 != -1)
      flag2 = true
    end

    #double move
    if flag2 == true
      if (abs(sourcex - targetx) == 1 || abs(sourcex - targetx) == 0) && ((abs(sourcey - targety) == 1) || (abs(sourcey2 - targety) == 0))
        if (abs(targetx2 - targetx) == 1 || abs(targetx2 - targetx) == 0) && ((abs(targety2 - targety) == 1) || (abs(targety2 - targety2) == 0))
          return true
        end
      end

    else
      if (abs(sourcex - targetx) == 2 || abs(soucey - targety) == 2)
        return true
      end
    end

  end #lionValidate end

  #case 16 sideMoverValidate found in case 10.2

  #case 16.2 freeBoarValidate #bishop and side mover
  function freeBoarValidate(team,sourcex,sourcey,targetx,targety)

    #check if moves like bishop
    if (moveDiagonalValidate(sourcex,sourcey,targetx,targety) == true)
      return true
    end

    #check if it moves like side mover
    if (sourcex == targetx)
      if (sourcey > targety) #moveLeftValidate
        if (moveLeftValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      else #moveRightValidate
        if (moveRightValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      end
    elseif (sourcey == targety)
      if (abs(sourcex - targetx) == 1)
        return true
      end
    end

    return false
  end #freeBoarValidate end

  #case 17 goBetweenValidate
  function goBetweenValidate(team,sourcex,sourcey,targetx,targety)
    if (sourcey == targety)
      if (abs(sourcey - targety) == 1)
        return true
      end
    end

    return false
  end #goBetweenValidate end

  #case 17.2 drunk elephant  implemented case 12

  #case 18 blindTigerValidate
  function blindTigerValidate(team,sourcex,sourcey,targetx,targety)

    if ((abs(sourcex - targetx) == 1) || (abs(sourcex - targetx) == 0))
      if ((abs(sourcey - targety) == 1) || (abs(sourcey - targety) == 0))
        if (team == 'w')
          if (sourcex != targetx - 1)
            return true #move down 1
          end
        elseif (team == 'b')
          if (sourcex != targetx + 1)
            return true #move up 1
          end
        end
      end
    end

    return false
  end #blindTigerValidate end

  #case 18.2 flygingstagvalidate
  function flyingStagValidate(team,sourcex,sourcey,targetx,targety)
    #moves like king
    if (moveAdjacentValidate(sourcex,sourcey,targetx,targety) == true)
      return true
    end

    #moves like vertical mover
    if (sourcey == targety)
      if ( sourcex < targetx) #move down
        if (moveDownValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      else #move up
        if (moveUpValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      end
    end

    return false
  end #flyingStagValidate end

  #case 19 queenValidate #moves like rook and bhishop
  function queenValidate(team,sourcex,sourcey,targetx,targety)

    #check if moves like bishop
    if (moveDiagonalValidate(sourcex,sourcey,targetx,targety) == true)
      return true
    end

    #moves like rook
    if (moveOrthogonalValidate(sourcex,sourcey,targetx,targety) == true)
      return true
    end

    return false
  end #queenValidate end

  #case 20 verticalMoverValidate
  function verticalMoverValidate(team,sourcex,sourcey,targetx,targety)

    #move up and down
    if (sourcey == targety)
      if (sourcex < targetx) #move down
        if (moveDownValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      else #move up
        if (moveUpValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      end
    end

    #left or right by 1
    if (sourcex == targetx)
      if (abs(sourcey - targety) == 1)
        return true
      end
    end

    return false
  end #verticalMoverValidate end

  #case 20.2 flyingOxValidate #moves like bishop and vertical mover
  function flyingOxValidate(team,sourcex,sourcey,targetx,targety)

    #check if moves like bishop
    if (moveDiagonalValidate(sourcex,sourcey,targetx,targety) == true)
      return true
    end

    #moves like vertical mover
    if (sourcey == targety)
      if (sourcex < targetx) #move down
        if (moveDownValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      else #move up
        if (moveUpValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      end
    end

    #left or right by 1
    if (sourcex == targetx)
      if (abs(sourcey - targety) == 1)
        return true
      end
    end

    return false
  end #flyingOxValidate

  #case 21 Phoenix
  function phoenixValidate(team,sourcex,sourcey,targetx,targety)

    if (abs(sourcex - targetx) == 2) && (abs(sourcey - targety) == 2) #jumps
      return true
    end

    if (sourcex == targetx) #up down left right by 1
      if (abs(sourcey - targety) == 1)
        return true
      end
    elseif (sourcey == targety)
      if (abs(sourcex - targetx) == 1)
        return true
      end
    end

    return false
  end #phoenixValidate end

  #case 21.2 queen implemented case 19

  #case 5 kirinValidate
  function kirinValidate(team,sourcex,sourcey,targetx,targety)

    if (abs(sourcex - targetx) == 1) && (abs(sourcey - targety) == 1) #jumps
      return true
    end

    if (sourcex == targetx) #left or right
      if (abs(sourcey - targety) == 2)
        return true
      end
    elseif (sourcey == targety) #up or down
      if (abs(sourcex - targety) == 2)
        return
      end
    end

    return false
  end #kirinValidate end

  #case 4.2 whiteHorseValidate
  function whiteHorseValidate(team,sourcex,sourcey,targetx,targety)

    if (sourcey == targety)
      #moving up
      if (sourcex > targetx)
        if (moveUpValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      #moving down
      else
        if (moveDownValidate(sourcex,sourcey,targetx,targety) == true)
          return true
        end
      end
    elseif (abs(sourcex - targetx) == abs(sourcey - targety)) #diagonal
      if (team == 'w') && (sourcex < targetx)
        if (sourcey < targety) #diagonalDownRight
          if (diagonalDownRightValidate(sourcex,sourcey,targetx,targety) == true)
            return true
          end
        else #diagonalDownLeft
          if (diagonalDownLeftValidate(sourcex,sourcey,targetx,targety) == true)
            return true
          end
        end
      elseif (team == 'b') && (sourcex > targetx)
        if (sourcey < targety) #diagonalUpRight
          if (diagonalUpRightValidate(sourcex,sourcey,targetx,targety) == true)
            return true
          end
        else #diagonalUpLeft
          if (diagonalUpLeftValidate(sourcex,sourcey,targetx,targety) == true)
            return true
          end
        end
      end
    end

    return false
  end #whiteHorseValidate end

end #validateMod

#<====== MISSING LIONVALIDATE AND LION MOVE FOR SOARING EAGLE =====> #

  maxMove = parse(get(SQLite.query(db, """SELECT max("move_number") from moves;""")[1,1]))
  validSoFar = true
  badMove=0
  metaMove = SQLite.query(db, """SELECT value FROM meta WHERE key = "type" """)
  for x in 1:maxMove #access each row of database
    dataMove = SQLite.query(db, """SELECT move_number, move_type, sourcex, sourcey, targetx, targety, targetx2, targety2, option, i_am_cheating FROM moves WHERE "move_number" = $x""" )
    if (!isnull(dataMove[1,5]) && !isnull(dataMove[1,3]) && !isnull(dataMove[1,6])) #targetx and targety not null

      gameType = get(metaMove[1,1])
      moveType = !isnull(dataMove[1,2]) ? get(dataMove[1,2]) : -1
      sourcex = !isnull(dataMove[1,3]) ? get(dataMove[1,3]) : -1
      sourcey = !isnull(dataMove[1,4]) ? get(dataMove[1,4]) : -1
      targetx = !isnull(dataMove[1,5]) ? get(dataMove[1,5]) : -1
      targety = !isnull(dataMove[1,6]) ? get(dataMove[1,6]) : -1
      targetx2 = !isnull(dataMove[1,7]) ? get(dataMove[1,7]) : -1
      targety2 = !isnull(dataMove[1,8]) ? get(dataMove[1,8]) : -1
      (unitType,team) = get(board,(sourcex,sourcey),('x','x'))
      if (moveValidate(unitType, moveType, gameType, team, sourcex, sourcey, targetx, targety, targetx2, targety2))
        #validSoFar

      else
        validSoFar = false
        badMove=get(dataMove[1,1])
        println(badMove)
      end
      board[(targetx,targety)] = deepcopy(board[(sourcex,sourcey)]) #Updates the board before next move
      delete!(board,(sourcex,sourcey))
    end
end
if !validSoFar #Ensures print only happens at end
  println(badMove)
else
  println(0)
end

max_stars_repo_name: kinoroy/Shogi

max_issues_repo_name: kinoroy/Shulia

max_forks_repo_name: kinoroy/Shulia


text:
a = [-0.1,0.2,0.3]
A = BandedMatrix(-2 => Vcat([1], Fill(1,∞)),  0 => Vcat(a, Fill(0,∞)), 1 => Vcat([1/4], Fill(1/4,∞)))
ℓ = λ -> ql(A-λ*I).L[1,1]
x,y = range(-2,2; length=200),range(-2,2;length=200)
    z = abs.(ℓ.(x' .+ y.*im))
    contour(x,y,z; nlevels=50, title="z^2 + 1/(4z)")
     scatter!(eigvals(Matrix(A[1:100,1:100])))
     scatter!(eigvals(Matrix(A[1:100,1:100]')))

a = [-0.5,0.2,0.3]
A = BandedMatrix(-2 => Vcat([1], Fill(1,∞)),  0 => Vcat(a, Fill(0,∞)), 1 => Vcat([1/4], Fill(1/4,∞)))
ℓ = λ -> ql(A-λ*I).L[1,1]
x,y = range(-2,2; length=200),range(-2,2;length=200)
    z = abs.(ℓ.(x' .+ y.*im))
    contour(x,y,z; nlevels=50, title="z^2 + 1/(4z)")
     scatter!(eigvals(Matrix(A[1:100,1:100])); label="A finite section")
     scatter!(eigvals(Matrix(A[1:100,1:100]')); label="A' finite section")

a = [-0.5,0.2,0.5]
A = BandedMatrix(-2 => Vcat([1], Fill(1,∞)),  0 => Vcat(a, Fill(0,∞)), 1 => Vcat([1/4], Fill(1/4,∞)))
ℓ = λ -> ql(A-λ*I).L[1,1]
x,y = range(-2,2; length=200),range(-2,2;length=200)
    z = abs.(ℓ.(x' .+ y.*im))
    contour(x,y,z; nlevels=50, title="z^2 + 1/(4z)")
    scatter!(eigvals(Matrix(A[1:100,1:100]')); label="A' finite section")
    scatter!(eigvals(Matrix(A[1:100,1:100])); label="A finite section")

x = range(-2,2,length=1_000)
    plot(x, abs.(ℓ.(x.+eps()im)); label="abs(L[1,1])")

a = [-0.1,0.2,0.3]
A = BandedMatrix(-2 => Vcat([1], Fill(1,∞)),  0 => Vcat(a, Fill(0,∞)), 1 => Vcat([1/4], Fill(1/4,∞)))
Q,L = ql(A-(0.5+0.1im)*I)
Q[1:10,1:11]*L[1:11,1:10] - (A-(0.5+0.1im)*I)[1:10,1:10]

ℓ(0.5+0.1im)
ql((A-(0.5+0.1im)*I)[1:1000,1:1000]).L[1,1]


text:
module euGENE

using Requests, DataFrames

export getGenBankIDs, getGenPeptIDs, getEntrezFasta, getUniProtFasta, parseFasta, translate, translateAll, parseORFs, reverseComplement, switchThyUra, calculateMW, cutAfter, cutBefore, calculateCharge, residueMass, bIons, yIons, fragmentationTable, fragmentationDataFrame, getResidues

RNA_CODONS = {"UUU" => 'F', "UUC" => 'F', "UUA" => 'L', "UUG" => 'L',
              "UCU" => 'S', "UCC" => 'S', "UCA" => 'S', "UCG" => 'S',
              "UAU" => 'Y', "UAC" => 'Y', "UAA" => '-', "UAG" => '-',
              "UGU" => 'C', "UGC" => 'C', "UGA" => '-', "UGG" => 'W',
              "CUU" => 'L', "CUC" => 'L', "CUA" => 'L', "CUG" => 'L',
              "CCU" => 'P', "CCC" => 'P', "CCA" => 'P', "CCG" => 'P',
              "CAU" => 'H', "CAC" => 'H', "CAA" => 'Q', "CAG" => 'Q',
              "CGU" => 'R', "CGC" => 'R', "CGA" => 'R', "CGG" => 'R',
              "AUU" => 'I', "AUC" => 'I', "AUA" => 'I', "AUG" => 'M',
              "ACU" => 'T', "ACC" => 'T', "ACA" => 'T', "ACG" => 'T',
              "AAU" => 'N', "AAC" => 'N', "AAA" => 'K', "AAG" => 'K',
              "AGU" => 'S', "AGC" => 'S', "AGA" => 'R', "AGG" => 'R',
              "GUU" => 'V', "GUC" => 'V', "GUA" => 'V', "GUG" => 'V',
              "GCU" => 'A', "GCC" => 'A', "GCA" => 'A', "GCG" => 'A',
              "GAU" => 'D', "GAC" => 'D', "GAA" => 'E', "GAG" => 'E',
              "GGU" => 'G', "GGC" => 'G', "GGA" => 'G', "GGG" => 'G'}

DNA_COMPLEMENT = {'A' => 'T', 'T' => 'A', 'C' => 'G', 'G' => 'C'}

# MW is for mass upon addition (water is already gone)
MW = {'A' => 71.03711, 'C' => 103.00919, 'D' => 115.02694, 'E' => 129.04259,
      'F' => 147.06841, 'G' => 57.02146, 'H' => 137.05891, 'I' => 113.08406,
      'K' => 128.09496, 'L' => 113.08406, 'M' => 131.04049, 'N' => 114.04293,
      'P' => 97.05276, 'Q' => 128.05858, 'R' => 156.10111, 'S' => 87.03203,
      'T' => 101.04768, 'V' => 99.06841, 'W' => 186.07931, 'Y' => 163.06333,
      # denotes carboxyamidomethylation, ~ denotes phosphorylation, @ denotes acetylation
      '#' => 57.02146, '~' => 79.96633, '@' => 42.01056}

ACID_CHAINS = {'C' => 8.14, 'D' => 3.71, 'E' => 4.15}
BASE_CHAINS = {'H' => 6.04, 'K' => 10.67, 'R' => 12.10}
SIDE_CHAINS = {'C' => 8.14, 'D' => 3.71, 'E' => 4.15, 'H' => 6.04, 'K' => 10.67, 'R' => 12.10}
AMINO_GROUPS = {'A' => 9.71, 'C' => 10.28, 'D' => 9.66, 'E' => 9.58,
                'F' => 9.09, 'G' => 9.58, 'H' => 9.09, 'I' => 9.60,
                'K' => 9.16, 'L' => 9.58, 'M' => 9.08, 'N' => 8.76,
                'P' => 10.47, 'Q' => 9.00, 'R' => 9.00, 'S' => 9.05,
                'T' => 8.96, 'V' => 9.52, 'W' => 9.34, 'Y' => 9.04}
CARBOXYL_GROUPS = {'A' => 2.33, 'C' => 1.91, 'D' => 1.95, 'E' => 2.16,
                   'F' => 2.18, 'G' => 2.34, 'H' => 1.70, 'I' => 2.26,
                   'K' => 2.15, 'L' => 2.32, 'M' => 2.16, 'N' => 2.16,
                   'P' => 1.95, 'Q' => 2.18, 'R' => 2.03, 'S' => 2.13,
                   'T' => 2.20, 'V' => 2.27, 'W' => 2.38, 'Y' => 2.24}

function getGenBankIDs(accession)
    return [n[5:end-5] for n in filter(x -> contains(x,"<Id>"),
    split(post("http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi";
               query = {"db" => "nucleotide", "term" => accession, "sort" => "relevance"}).data,r"\s"))]
end

function getGenPeptIDs(accession)
    return [n[5:end-5] for n in filter(x -> contains(x,"<Id>"),
    split(post("http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi";
               query = {"db" => "protein", "term" => accession, "sort" => "relevance"}).data,r"\s"))]
end

function getEntrezFasta(accession)
    return parseFasta(post("http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi";
                           query = {"db" => "protein", "id" => accession, "retmode" => "text", "rettype" => "fasta"}).data)
end

function getUniProtFasta(accession)
    return parseFasta(post("http://www.uniprot.org/uniprot/" * accession * ".fasta").data)
end

function parseFasta(fasta)
    return join(filter(x -> !contains(x,">"),split(fasta,"\n")))
end

# list comprehension version, uses same algorithm
function translate(rna)
    return join([RNA_CODONS[rna[n-2:n]] for n in 3:3:length(rna)])
end

# returns array of 5'->3' frame 1, 5'->3' frame 2, 5'->3' frame 3, 3'->5' frame 1, 3'->5' frame 2, 3'->5' frame 3
function translateAll(rna)
    return [translate(rna), translate(rna[2:end]), translate(rna[3:end]),
            translate(reverseComplement(rna)), translate(reverseComplement(rna)[2:end]), translate(reverseComplement(rna)[3:end])]
end

# list comprehension version, uses same algorithm
function parseORFs(translation)
    return [frame[search(frame,r"M.*")] for frame in filter(x -> contains(x,"M"),split(translation,"-"))]
end

function reverseComplement(sequence)
    if contains(sequence,"T")
        return join([DNA_COMPLEMENT[nucleotide] for nucleotide in reverse(sequence)])
    else
        return switchThyUra(join([DNA_COMPLEMENT[nucleotide] for nucleotide in reverse(switchThyUra(sequence))]))
    end
end

# if both T and U are in the sequence, it will still switch them (no judgments here)
function switchThyUra(sequence)
    return replace(replace(replace(uppercase(sequence),"T","~"),"U","T"),"~","U")
end

function calculateMW(sequence)
    return 18.053 + sum([residueMass(aa) for aa in collect(sequence)])
end

function calculatePI(sequence)
    pKas = [SIDE_CHAINS[n] for n in filter(x -> in(x,"CDEHKR"), collect(sequence))]
    push!(pKas,AMINO_GROUPS[sequence[1]])
    push!(pKas,CARBOXYL_GROUPS[sequence[end]])
    sort!(pKas)
    positive = 1 + count(x -> in(x,"HKR"), collect(sequence))
    return (pKas[positive] + pKas[positive+1]) / 2
end

function cutAfter(sequence,residues)
    fragments = [convert(ASCIIString,n) for n in matchall(Regex("[^$residues]*[$residues]"),sequence)]
    if length(join(fragments)) != length(sequence)
        push!(fragments,sequence[length(join(fragments))+1:end])
    end
    return fragments
end

function cutBefore(sequence,residues)
    fragments = [convert(ASCIIString,n) for n in matchall(Regex("[$residues][^$residues]*"),sequence)]
    if length(join(fragments)) != length(sequence)
        unshift!(fragments,sequence[1:length(sequence) - length(join(fragments))])
    end
    return fragments
end

function calculateCharge(sequence)
    return 1 + count(c -> in(c,"RK"), collect(sequence))
end

function residueMass(residue)
    return sum([MW[residue[n]] for n in 1:length(residue)])
end

function bIonMZ(sequence)
    return (sum([residueMass(aa) for aa in collect(sequence)]) + calculateCharge(sequence) * 1.00794) / calculateCharge(sequence)
end

function bIonMZ(sequence,charge)
    return (sum([residueMass(aa) for aa in collect(sequence)]) + charge * 1.00794) / charge
end

function yIonMZ(sequence)
    return (sum([residueMass(aa) for aa in collect(sequence)]) + calculateCharge(sequence) * 1.00794 + 18.01) / calculateCharge(sequence)
end

function yIonMZ(sequence,charge)
    return (sum([residueMass(aa) for aa in collect(sequence)]) + charge * 1.00794 + 18.01) / charge
end

function bIons(sequence)
    residues = getResidues(sequence)
    return [bIonMZ(stringJoin(residues[1:n])) for n in 1:1:length(residues)]
end

function bIons(sequence,charge)
    residues = getResidues(sequence)
    return [bIonMZ(stringJoin(residues[1:n]),charge) for n in 1:1:length(residues)]
end

function yIons(sequence)
    residues = getResidues(sequence)
    return [yIonMZ(stringJoin(reverse(residues[1:n]))) for n in 1:1:length(residues)]
end

function yIons(sequence,charge)
    residues = getResidues(sequence)
    return [yIonMZ(stringJoin(reverse(residues[1:n])),charge) for n in 1:1:length(residues)]
end

function fragmentationTable(sequence,charge)
    residues = getResidues(sequence)
    b = bIons(sequence,charge)
    y = reverse(yIons(sequence,charge))
    return [[residues[n],b[n],y[n]] for n in 1:1:length(sequence)]
end

function fragmentationDataFrame(sequence,charge)
    return DataFrame(Seq = getResidues(sequence), b = bIons(sequence,charge), y = reverse(yIons(sequence,charge)))
end

function getResidues(sequence)
    indices = Int[]
    for n in 1:length(sequence)
        if !in(sequence[n],"@~#")
            push!(indices,n)
        end
    end
    residues = [sequence[indices[n]:(indices[n+1]-1)] for n in 1:(length(indices)-1)]
    push!(residues,sequence[indices[end]:end])
    return residues
end

function stringJoin(array)
    buffer = ""
    for n in 1:length(array)
        buffer *= array[n]
    end
    return buffer
end

# module end
end


text:
using LightGraphs
using Distributions
using MonteCarloX
using Random
using StatsBase

# this runs 5s for run(1000,0.1,1.0,1e-3,1e5,1/1e-3,1000) as does my C code... :P
#
# TODO: reduce allocations? pass arrays?

function run(N::Int, p::Float64, m::Float64, h::Float64, T::Float64, T_therm::Float64, seed::Int; flag_fast = true)::Tuple{Float64,Float64,Int}
    rng = MersenneTwister(seed);

    # dimensionless rates
    mu     = 1.0
    lambda = mu * m

    system = construct_system(N, p, mu, lambda, h, seed, initial = "empty")

    T_total = T + T_therm;
    time = dt = dtime = 0.0;
    dtime_step = 10.0;
    N_active = N_meas = 0;
    avg_activity = avg_activity2 = 0.0;

    # TODO: think about abstract Gillespie: pass only system, rng, update(), external events as list
    # Q?? Can we do something similar for Metropolis?
    events_since_last_sum = 0
    while time < T_total
    # find next event (dt and time) to be updated
        dt, n = next(system.rates, system.sum_rates, rng)
        events_since_last_sum += 1
        if events_since_last_sum > 1000
            system.sum_rates = sum(system.rates)
            events_since_last_sum = 0
        end

        # measure observables in discrete time
        # use values from last step because they were valid inbetween
        dtime += dt;
        while dtime > dtime_step
            dtime     -= dtime_step
          # TODO: is this correct?
            if time > T_therm
                avg_activity  += N_active
                avg_activity2 += N_active * N_active
                N_meas        += 1
            end
        end
        # TODO: what is correct? T=1e4 with dt=10 -> N=1000 or 1001? here 1001

        # only after measurement evolve real time
        time += dt;

        N_active += update(system, n)
    end

    return avg_activity / N_meas, avg_activity2 / N_meas, N_meas
end

# crucial here is type-stability achieved with F1 and F2!!!
mutable struct ContactProcess{F1,F2}
    network::SimpleDiGraph{Int64}
    incoming_neighbors::F1
    outgoing_neighbors::F2
    neurons::Vector{Int8}
    active_incoming_neighbors::Vector{Int8}
    rates::ProbabilityWeights{Float64,Float64,Vector{Float64}}  # static list of rates
    mu::Float64
    lambda::Float64
    h::Float64
end

# here one could select initial condition empty, full, equilibrium?
function construct_system(N, p, mu, lambda, h, seed; initial = "empty")::ContactProcess
    # directed ER graph
    network = LightGraphs.SimpleGraphs.erdos_renyi(N, p, is_directed = true, seed = seed)
    incoming_neighbors(n) = inneighbors(network, n)
    outgoing_neighbors(n) = outneighbors(network, n)

    # initial conditions
    neurons = []
    active_incoming_neighbors = []
    if initial == "empty"
        neurons = zeros(Int8, N)
        active_incoming_neighbors = zeros(Int8, N);
        rates = ProbabilityWeights(ones(Float64, N) .* h);
    end
    if initial == "full"
        neurons = ones(Int8, N);
        active_incoming_neighbors = ones(Int8, N);
        rates = ProbabilityWeights(ones(Float64, N) .* mu);
    end

    return ContactProcess(network, incoming_neighbors, outgoing_neighbors, neurons, active_incoming_neighbors, rates, mu, lambda, h)
end

# system has to include list_rates and sum_ratesand neuon_active_neighbors, graph, indegree function etc
function update(system::ContactProcess, n::Int)::Int
    rate_n_old = system.rates[n];
    dstate     = 0;
    if system.neurons[n] == 0
        system.neurons[n] = 1;
        system.rates[n]   = system.mu;
        dstate            = 1;
    else
        system.neurons[n] = 0;
        system.rates[n]   = rate_activate(system, n)
        dstate            = -1;
    end

    # update rate list
    system.sum_rates += (system.rates[n] - rate_n_old);
    for nn in system.outgoing_neighbors(n)
        system.active_incoming_neighbors[nn] += dstate;
        if system.neurons[nn] == 0
            rate_nn_old       = system.rates[nn];
            system.rates[nn]  = rate_activate(system, nn)
            system.sum_rates += (system.rates[nn] - rate_nn_old);
        end
    end

    return dstate
end

@inline function rate_activate(system::ContactProcess, n::Int)::Float64
    return system.h + system.lambda * system.active_incoming_neighbors[n] / length(system.incoming_neighbors(n))
end

##############################run
# avg_activity, avg_activity2, T = run(Int(1e3),1.0,1e-3,1e-2,1e6,1e3,rng)
# println(avg_activity, " ", avg_activity2, " ", T)

# err = sqrt((avg_activity2-avg_activity*avg_activity)/(T-1))
# println(avg_activity, " ", err, " ", T)


text:
using RecipesBase

@recipe function f(lines::Vector{Line{T,2}}) where {T}
    xs = hcat(([line.p₁[1],line.p₂[1]] for line ∈ lines)...)
    ys = hcat(([line.p₁[2],line.p₂[2]] for line ∈ lines)...)
    linecolor --> :black
    legend := false
    xs,ys
end

@recipe function f(nwn::NWN_circuit)
    legend -> false
    grid -> false
    seriestype := :line
    @series begin
        linecolor := :black
        label := "Wires"
        nwn.lines
    end
    @series begin
        linecolor := :red
        label := "Electrode"
        linewidth --> 2
        nwn.elecs
    end
    @series begin
        linecolor := :blue
        label := "Ground"
        linewidth --> 2
        nwn.grnds
    end
end

@recipe function f(net::StickNetwork)
    legend -> false
    grid -> false
    seriestype := :line
    @series begin
        linecolor := :black
        net.lines
    end
end