From 86ec8eada3ea5209ca2072a639180a81aec8f0ec Mon Sep 17 00:00:00 2001 From: KINDNICK <68893236+KINDNICK@users.noreply.github.com> Date: Thu, 15 May 2025 00:31:20 +0900 Subject: [PATCH] =?UTF-8?q?ADD=20:=20=ED=95=B8=EB=93=9C=20=ED=8A=B8?= =?UTF-8?q?=EB=9E=98=ED=82=B9=20=ED=85=8C=EC=8A=A4=ED=8A=B8=20=EC=8A=A4?= =?UTF-8?q?=ED=81=AC=EB=A6=BD=ED=8A=B8=20=EC=97=85=EB=8D=B0=EC=9D=B4?= =?UTF-8?q?=ED=8A=B8?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Assets/External/BioIK.meta | 8 + Assets/External/BioIK/BioIK.cs | 299 + Assets/External/BioIK/BioIK.cs.meta | 19 + Assets/External/BioIK/Demo.meta | 9 + Assets/External/BioIK/Demo/Demo.unity | 3 + Assets/External/BioIK/Demo/Demo.unity.meta | 15 + .../External/BioIK/Demo/DemoSettings.lighting | 63 + .../BioIK/Demo/DemoSettings.lighting.meta | 8 + Assets/External/BioIK/Demo/Materials.meta | 9 + Assets/External/BioIK/Demo/Materials/Blue.mat | 139 + .../BioIK/Demo/Materials/Blue.mat.meta | 15 + Assets/External/BioIK/Demo/Materials/Box.mat | 139 + .../BioIK/Demo/Materials/Box.mat.meta | 15 + .../External/BioIK/Demo/Materials/Green.mat | 138 + .../BioIK/Demo/Materials/Green.mat.meta | 16 + .../External/BioIK/Demo/Materials/Ground.mat | 139 + .../BioIK/Demo/Materials/Ground.mat.meta | 15 + Assets/External/BioIK/Demo/Materials/Red.mat | 139 + .../BioIK/Demo/Materials/Red.mat.meta | 15 + .../External/BioIK/Demo/Materials/Target.mat | 139 + .../BioIK/Demo/Materials/Target.mat.meta | 15 + Assets/External/BioIK/Demo/Materials/Wall.mat | 139 + .../BioIK/Demo/Materials/Wall.mat.meta | 15 + Assets/External/BioIK/Demo/Meshes.meta | 9 + Assets/External/BioIK/Demo/Meshes/Cone.asset | 3 + .../BioIK/Demo/Meshes/Cone.asset.meta | 15 + Assets/External/BioIK/Demo/Prefabs.meta | 9 + .../External/BioIK/Demo/Prefabs/Target.prefab | 3 + .../BioIK/Demo/Prefabs/Target.prefab.meta | 15 + Assets/External/BioIK/Demo/Robot Kyle.meta | 9 + .../Demo/Robot Kyle/Controller.controller | 867 + .../Robot Kyle/Controller.controller.meta | 13 + .../BioIK/Demo/Robot Kyle/Materials.meta | 9 + .../Demo/Robot Kyle/Materials/Robot_Color.mat | 140 + .../Robot Kyle/Materials/Robot_Color.mat.meta | 13 + .../External/BioIK/Demo/Robot Kyle/Model.meta | 9 + .../Demo/Robot Kyle/Model/Robot Kyle.fbx | 3 + .../Demo/Robot Kyle/Model/Robot Kyle.fbx.meta | 829 + .../External/BioIK/Demo/Robot Kyle/Run.anim | 51434 ++++++++++++++++ .../BioIK/Demo/Robot Kyle/Run.anim.meta | 13 + .../BioIK/Demo/Robot Kyle/Textures.meta | 9 + .../Demo/Robot Kyle/Textures/Robot_Color.tga | 3 + .../Robot Kyle/Textures/Robot_Color.tga.meta | 64 + .../Demo/Robot Kyle/Textures/Robot_Normal.tga | 3 + .../Robot Kyle/Textures/Robot_Normal.tga.meta | 64 + Assets/External/BioIK/Demo/Scripts.meta | 9 + .../BioIK/Demo/Scripts/CameraController.cs | 83 + .../Demo/Scripts/CameraController.cs.meta | 19 + Assets/External/BioIK/Demo/Scripts/Demo.cs | 84 + .../External/BioIK/Demo/Scripts/Demo.cs.meta | 19 + .../External/BioIK/Demo/Scripts/MouseDrag.cs | 62 + .../BioIK/Demo/Scripts/MouseDrag.cs.meta | 19 + Assets/External/BioIK/Demo/Scripts/Rotate.cs | 13 + .../BioIK/Demo/Scripts/Rotate.cs.meta | 19 + Assets/External/BioIK/Demo/Textures.meta | 9 + Assets/External/BioIK/Demo/Textures/Steel.jpg | 3 + .../BioIK/Demo/Textures/Steel.jpg.meta | 66 + Assets/External/BioIK/Editor.meta | 9 + Assets/External/BioIK/Editor/BioIKEditor.cs | 840 + .../External/BioIK/Editor/BioIKEditor.cs.meta | 19 + Assets/External/BioIK/Helpers.meta | 9 + Assets/External/BioIK/Helpers/Enums.cs | 5 + Assets/External/BioIK/Helpers/Enums.cs.meta | 19 + Assets/External/BioIK/Helpers/Evolution.cs | 6105 ++ .../External/BioIK/Helpers/Evolution.cs.meta | 19 + Assets/External/BioIK/Helpers/Model.cs | 577 + Assets/External/BioIK/Helpers/Model.cs.meta | 19 + Assets/External/BioIK/Helpers/Utility.cs | 177 + Assets/External/BioIK/Helpers/Utility.cs.meta | 19 + Assets/External/BioIK/ReadMe | 95 + Assets/External/BioIK/ReadMe.meta | 15 + Assets/External/BioIK/Setup.meta | 9 + Assets/External/BioIK/Setup/BioJoint.cs | 538 + Assets/External/BioIK/Setup/BioJoint.cs.meta | 19 + Assets/External/BioIK/Setup/BioObjective.cs | 82 + .../External/BioIK/Setup/BioObjective.cs.meta | 19 + Assets/External/BioIK/Setup/BioSegment.cs | 75 + .../External/BioIK/Setup/BioSegment.cs.meta | 19 + Assets/External/BioIK/Setup/Objectives.meta | 9 + .../BioIK/Setup/Objectives/Displacement.cs | 47 + .../Setup/Objectives/Displacement.cs.meta | 19 + .../BioIK/Setup/Objectives/Distance.cs | 133 + .../BioIK/Setup/Objectives/Distance.cs.meta | 19 + .../BioIK/Setup/Objectives/JointValue.cs | 131 + .../BioIK/Setup/Objectives/JointValue.cs.meta | 19 + .../External/BioIK/Setup/Objectives/LookAt.cs | 123 + .../BioIK/Setup/Objectives/LookAt.cs.meta | 19 + .../BioIK/Setup/Objectives/Orientation.cs | 106 + .../Setup/Objectives/Orientation.cs.meta | 19 + .../BioIK/Setup/Objectives/Position.cs | 99 + .../BioIK/Setup/Objectives/Position.cs.meta | 19 + .../BioIK/Setup/Objectives/Projection.cs | 231 + .../BioIK/Setup/Objectives/Projection.cs.meta | 19 + Assets/Scripts/HandTracking.meta | 8 + Assets/Scripts/HandTracking/HandTracking.meta | 8 + .../Scripts/HandTracking/HandTracking.unity | 3 + .../HandTracking/HandTracking.unity.meta | 7 + .../HandTracking/GlobalVolumeProfile.asset | 3 + .../GlobalVolumeProfile.asset.meta | 8 + Assets/Scripts/HandTracking/PuppetHand.fbx | 3 + .../Scripts/HandTracking/PuppetHand.fbx.meta | 107 + 101 files changed, 65299 insertions(+) create mode 100644 Assets/External/BioIK.meta create mode 100644 Assets/External/BioIK/BioIK.cs create mode 100644 Assets/External/BioIK/BioIK.cs.meta create mode 100644 Assets/External/BioIK/Demo.meta create mode 100644 Assets/External/BioIK/Demo/Demo.unity create mode 100644 Assets/External/BioIK/Demo/Demo.unity.meta create mode 100644 Assets/External/BioIK/Demo/DemoSettings.lighting create mode 100644 Assets/External/BioIK/Demo/DemoSettings.lighting.meta create mode 100644 Assets/External/BioIK/Demo/Materials.meta create mode 100644 Assets/External/BioIK/Demo/Materials/Blue.mat create mode 100644 Assets/External/BioIK/Demo/Materials/Blue.mat.meta create mode 100644 Assets/External/BioIK/Demo/Materials/Box.mat create mode 100644 Assets/External/BioIK/Demo/Materials/Box.mat.meta create mode 100644 Assets/External/BioIK/Demo/Materials/Green.mat create mode 100644 Assets/External/BioIK/Demo/Materials/Green.mat.meta create mode 100644 Assets/External/BioIK/Demo/Materials/Ground.mat create mode 100644 Assets/External/BioIK/Demo/Materials/Ground.mat.meta create mode 100644 Assets/External/BioIK/Demo/Materials/Red.mat create mode 100644 Assets/External/BioIK/Demo/Materials/Red.mat.meta create mode 100644 Assets/External/BioIK/Demo/Materials/Target.mat create mode 100644 Assets/External/BioIK/Demo/Materials/Target.mat.meta create mode 100644 Assets/External/BioIK/Demo/Materials/Wall.mat create mode 100644 Assets/External/BioIK/Demo/Materials/Wall.mat.meta create mode 100644 Assets/External/BioIK/Demo/Meshes.meta create mode 100644 Assets/External/BioIK/Demo/Meshes/Cone.asset create mode 100644 Assets/External/BioIK/Demo/Meshes/Cone.asset.meta create mode 100644 Assets/External/BioIK/Demo/Prefabs.meta create mode 100644 Assets/External/BioIK/Demo/Prefabs/Target.prefab create mode 100644 Assets/External/BioIK/Demo/Prefabs/Target.prefab.meta create mode 100644 Assets/External/BioIK/Demo/Robot Kyle.meta create mode 100644 Assets/External/BioIK/Demo/Robot Kyle/Controller.controller create mode 100644 Assets/External/BioIK/Demo/Robot Kyle/Controller.controller.meta create mode 100644 Assets/External/BioIK/Demo/Robot Kyle/Materials.meta create mode 100644 Assets/External/BioIK/Demo/Robot Kyle/Materials/Robot_Color.mat create mode 100644 Assets/External/BioIK/Demo/Robot Kyle/Materials/Robot_Color.mat.meta create mode 100644 Assets/External/BioIK/Demo/Robot Kyle/Model.meta create mode 100644 Assets/External/BioIK/Demo/Robot Kyle/Model/Robot Kyle.fbx create mode 100644 Assets/External/BioIK/Demo/Robot Kyle/Model/Robot Kyle.fbx.meta create mode 100644 Assets/External/BioIK/Demo/Robot Kyle/Run.anim create mode 100644 Assets/External/BioIK/Demo/Robot Kyle/Run.anim.meta create mode 100644 Assets/External/BioIK/Demo/Robot Kyle/Textures.meta create mode 100644 Assets/External/BioIK/Demo/Robot Kyle/Textures/Robot_Color.tga create mode 100644 Assets/External/BioIK/Demo/Robot Kyle/Textures/Robot_Color.tga.meta create mode 100644 Assets/External/BioIK/Demo/Robot Kyle/Textures/Robot_Normal.tga create mode 100644 Assets/External/BioIK/Demo/Robot Kyle/Textures/Robot_Normal.tga.meta create mode 100644 Assets/External/BioIK/Demo/Scripts.meta create mode 100644 Assets/External/BioIK/Demo/Scripts/CameraController.cs create mode 100644 Assets/External/BioIK/Demo/Scripts/CameraController.cs.meta create mode 100644 Assets/External/BioIK/Demo/Scripts/Demo.cs create mode 100644 Assets/External/BioIK/Demo/Scripts/Demo.cs.meta create mode 100644 Assets/External/BioIK/Demo/Scripts/MouseDrag.cs create mode 100644 Assets/External/BioIK/Demo/Scripts/MouseDrag.cs.meta create mode 100644 Assets/External/BioIK/Demo/Scripts/Rotate.cs create mode 100644 Assets/External/BioIK/Demo/Scripts/Rotate.cs.meta create mode 100644 Assets/External/BioIK/Demo/Textures.meta create mode 100644 Assets/External/BioIK/Demo/Textures/Steel.jpg create mode 100644 Assets/External/BioIK/Demo/Textures/Steel.jpg.meta create mode 100644 Assets/External/BioIK/Editor.meta create mode 100644 Assets/External/BioIK/Editor/BioIKEditor.cs create mode 100644 Assets/External/BioIK/Editor/BioIKEditor.cs.meta create mode 100644 Assets/External/BioIK/Helpers.meta create mode 100644 Assets/External/BioIK/Helpers/Enums.cs create mode 100644 Assets/External/BioIK/Helpers/Enums.cs.meta create mode 100644 Assets/External/BioIK/Helpers/Evolution.cs create mode 100644 Assets/External/BioIK/Helpers/Evolution.cs.meta create mode 100644 Assets/External/BioIK/Helpers/Model.cs create mode 100644 Assets/External/BioIK/Helpers/Model.cs.meta create mode 100644 Assets/External/BioIK/Helpers/Utility.cs create mode 100644 Assets/External/BioIK/Helpers/Utility.cs.meta create mode 100644 Assets/External/BioIK/ReadMe create mode 100644 Assets/External/BioIK/ReadMe.meta create mode 100644 Assets/External/BioIK/Setup.meta create mode 100644 Assets/External/BioIK/Setup/BioJoint.cs create mode 100644 Assets/External/BioIK/Setup/BioJoint.cs.meta create mode 100644 Assets/External/BioIK/Setup/BioObjective.cs create mode 100644 Assets/External/BioIK/Setup/BioObjective.cs.meta create mode 100644 Assets/External/BioIK/Setup/BioSegment.cs create mode 100644 Assets/External/BioIK/Setup/BioSegment.cs.meta create mode 100644 Assets/External/BioIK/Setup/Objectives.meta create mode 100644 Assets/External/BioIK/Setup/Objectives/Displacement.cs create mode 100644 Assets/External/BioIK/Setup/Objectives/Displacement.cs.meta create mode 100644 Assets/External/BioIK/Setup/Objectives/Distance.cs create mode 100644 Assets/External/BioIK/Setup/Objectives/Distance.cs.meta create mode 100644 Assets/External/BioIK/Setup/Objectives/JointValue.cs create mode 100644 Assets/External/BioIK/Setup/Objectives/JointValue.cs.meta create mode 100644 Assets/External/BioIK/Setup/Objectives/LookAt.cs create mode 100644 Assets/External/BioIK/Setup/Objectives/LookAt.cs.meta create mode 100644 Assets/External/BioIK/Setup/Objectives/Orientation.cs create mode 100644 Assets/External/BioIK/Setup/Objectives/Orientation.cs.meta create mode 100644 Assets/External/BioIK/Setup/Objectives/Position.cs create mode 100644 Assets/External/BioIK/Setup/Objectives/Position.cs.meta create mode 100644 Assets/External/BioIK/Setup/Objectives/Projection.cs create mode 100644 Assets/External/BioIK/Setup/Objectives/Projection.cs.meta create mode 100644 Assets/Scripts/HandTracking.meta create mode 100644 Assets/Scripts/HandTracking/HandTracking.meta create mode 100644 Assets/Scripts/HandTracking/HandTracking.unity create mode 100644 Assets/Scripts/HandTracking/HandTracking.unity.meta create mode 100644 Assets/Scripts/HandTracking/HandTracking/GlobalVolumeProfile.asset create mode 100644 Assets/Scripts/HandTracking/HandTracking/GlobalVolumeProfile.asset.meta create mode 100644 Assets/Scripts/HandTracking/PuppetHand.fbx create mode 100644 Assets/Scripts/HandTracking/PuppetHand.fbx.meta diff --git a/Assets/External/BioIK.meta b/Assets/External/BioIK.meta new file mode 100644 index 00000000..de379259 --- /dev/null +++ b/Assets/External/BioIK.meta @@ -0,0 +1,8 @@ +fileFormatVersion: 2 +guid: 5922c52739f8e0346b7b85d9c293f538 +folderAsset: yes +DefaultImporter: + externalObjects: {} + userData: + assetBundleName: + assetBundleVariant: diff --git a/Assets/External/BioIK/BioIK.cs b/Assets/External/BioIK/BioIK.cs new file mode 100644 index 00000000..052bdd3f --- /dev/null +++ b/Assets/External/BioIK/BioIK.cs @@ -0,0 +1,299 @@ +using UnityEngine; +using System.Collections.Generic; + +namespace BioIK { + + //[ExecuteInEditMode] + [DisallowMultipleComponent] + public class BioIK : MonoBehaviour { + + //public bool SolveInEditMode = false; + + [SerializeField] private bool UseThreading = true; + + [SerializeField] private int Generations = 2; + [SerializeField] private int PopulationSize = 50; + [SerializeField] private int Elites = 2; + + public float Smoothing = 0.5f; + public float AnimationWeight = 0f; + public float AnimationBlend = 0f; + public MotionType MotionType = MotionType.Instantaneous; + public float MaximumVelocity = 3f; + public float MaximumAcceleration = 3f; + + public List Segments = new List(); + + public BioSegment Root = null; + public Evolution Evolution = null; + public double[] Solution = null; + + private bool Destroyed = false; + + //Custom Inspector Helpers + public BioSegment SelectedSegment = null; + public Vector2 Scroll = Vector2.zero; + + void Awake() { + Refresh(); + } + + void Start() { + + } + + void OnDestroy() { + Destroyed = true; + DeInitialise(); + Utility.Cleanup(transform); + } + + void OnEnable() { + Initialise(); + } + + void OnDisable() { + DeInitialise(); + } + + private void Initialise() { + if(Evolution == null) { + Evolution = new Evolution(new Model(this), PopulationSize, Elites, UseThreading); + } + } + + private void DeInitialise() { + if(Evolution != null) { + Evolution.Kill(); + Evolution = null; + } + } + + void Update() { + PrecaptureAnimation(Root); + } + + void LateUpdate() { + PostcaptureAnimation(Root); + + UpdateData(Root); + + for(int i=0; i GetChain(Transform start, Transform end) { + BioSegment a = FindSegment(start); + BioSegment b = FindSegment(end); + if(a == null || b == null) { + Debug.Log("Could not generate chain for given transforms"); + return null; + } + return GetChain(a, b); + } + + public List GetChain(BioSegment start, BioSegment end) { + List chain = new List(); + BioSegment segment = end; + while(true) { + chain.Add(segment); + if(segment.Transform == transform || segment.Parent == null) { + break; + } else { + segment = segment.Parent; + } + } + chain.Reverse(); + return chain; + } + + public void UpdateData(BioSegment segment) { + if(segment.Joint != null) { + if(segment.Joint.enabled) { + segment.Joint.UpdateData(); + } + } + for(int i=0; i index) { + if(ModelGO != null) { + ModelGO.SetActive(false); + } + ModelCharacter = Models[index]; + ModelGO = ModelCharacter.transform.root.gameObject; + ModelGO.SetActive(true); + UpdateMotionType(); + } + } + + public bool Ignore(GameObject go) { + for(int i=0; i()) { + t.GetComponent().hideFlags = HideFlags.None; + } + foreach(BioObjective o in t.GetComponents()) { + o.hideFlags = HideFlags.None; + } + if(t.GetComponent()) { + t.GetComponent().hideFlags = HideFlags.None; + } + for(int i=0; i()) { + t.GetComponent().hideFlags = HideFlags.HideInInspector; + } + foreach(BioObjective o in t.GetComponents()) { + o.hideFlags = HideFlags.HideInInspector; + } + if(t.GetComponent()) { + t.GetComponent().hideFlags = HideFlags.HideInInspector; + } + for(int i=0; i 0) { + SetGUIColor(Color13); + using(new EditorGUILayout.VerticalScope ("Box")) { + int width = 10*(indent-1); + EditorGUILayout.LabelField("", GUILayout.Width(width)); + } + } + } + + GUI.skin.button.alignment = TextAnchor.MiddleLeft; + if(segment == Target.SelectedSegment) { + SetGUIColor(Color5); + } else { + SetGUIColor(Color.Lerp(Color4, Color8, (float)indent / (float)maxIndent)); + } + if(GUILayout.Button(segment.Transform.name, GUILayout.Height(25f), GUILayout.ExpandWidth(true))) { + if(Target.SelectedSegment == segment) { + Target.SelectedSegment = null; + ChosingObjectiveType = false; + } else { + Target.SelectedSegment = segment; + } + } + + EditorGUILayout.EndHorizontal(); + + if(Target.SelectedSegment == segment) { + InspectSegment(segment); + } else { + GUILayout.BeginHorizontal(); + GUILayout.FlexibleSpace(); + if(segment.Joint != null) { + SetGUIColor(Color6); + GUILayout.Box(" Joint "); + } + foreach(BioObjective objective in segment.Objectives) { + SetGUIColor(Color9); + GUILayout.Box(" " + objective.GetObjectiveType().ToString() + " "); + } + GUILayout.FlexibleSpace(); + GUILayout.EndHorizontal(); + } + } + + for(int i=0; i Pool = new List(); //Selection pool for recombination + private int PoolCount; //Current size of the selection pool + private double[] Probabilities; //Current probabilities for selection + private double Gene; //Simple storage variable #1 + private double Weight; //Simple storage variable #2 + private bool[] Constrained; + + //Variables for elitism exploitation + private bool UseThreading; + private bool Evolving = false; + private bool[] Improved = null; + private Model[] Models = null; + private BFGS[] Optimisers = null; + + //Threading + private bool ThreadsRunning = false; + private ManualResetEvent[] Handles = null; + private Thread[] Threads = null; + private bool[] Work = null; + + //Variables for optimisation queries + private double[] Solution; //Evolutionary solution + private double Fitness; //Evolutionary fitness + + private bool Killed = false; + + //Initialises the algorithm + public Evolution(Model model, int populationSize, int elites, bool useThreading) { + Model = model; + PopulationSize = populationSize; + Elites = elites; + Dimensionality = model.GetDoF(); + UseThreading = useThreading; + + Population = new Individual[PopulationSize]; + Offspring = new Individual[PopulationSize]; + for(int i=0; i Models[index].ComputeLoss(x), y => Models[index].ComputeGradient(y, 1e-5)); + } + + if(UseThreading) { + //Start Threads + ThreadsRunning = true; + Work = new bool[Elites]; + Handles = new ManualResetEvent[Elites]; + Threads = new Thread[Elites]; + for(int i=0; i SurviveThread(index)); + //Threads[index].Start(); + } + } else { + ThreadsRunning = false; + } + } + + ~Evolution() { + Kill(); + } + + public void Kill() { + if(Killed) { + return; + } + Killed = true; + if(UseThreading) { + //Stop Threads + ThreadsRunning = false; + for(int i=0; i 0) { + Individual parentA = Select(Pool); + Individual parentB = Select(Pool); + Individual prototype = Select(Pool); + + //Recombination and Adoption + Reproduce(Offspring[i], parentA, parentB, prototype); + + //Pre-Selection Niching + if(Offspring[i].Fitness < parentA.Fitness) { + Pool.Remove(parentA); + PoolCount -= 1; + } + if(Offspring[i].Fitness < parentB.Fitness) { + Pool.Remove(parentB); + PoolCount -= 1; + } + } else { + //Fill the population + Reroll(Offspring[i]); + } + } + + //Finish evolving + Evolving = false; + + //Wait for threads to finish + //System.DateTime timer = Utility.GetTimestamp(); + while(HasWork()) { + //Wait one cycle + /* + if(Utility.GetElapsedTime(timer) > 0.5f) { + Debug.Log("!!! !!! !!! FIXING CRASH !!! !!! !!!"); + for(int i=0; i 0) { + Individual parentA = Select(Pool); + Individual parentB = Select(Pool); + Individual prototype = Select(Pool); + + //Recombination and Adoption + Reproduce(Offspring[i], parentA, parentB, prototype); + + //Pre-Selection Niching + if(Offspring[i].Fitness < parentA.Fitness) { + Pool.Remove(parentA); + PoolCount -= 1; + } + if(Offspring[i].Fitness < parentB.Fitness) { + Pool.Remove(parentB); + PoolCount -= 1; + } + } else { + //Fill the population + Reroll(Offspring[i]); + } + } + double duration = Utility.GetElapsedTime(timestamp); + + //Exploit elites sequentially + for(int i=0; i UpperBounds[i]) { + offspring.Genes[i] = UpperBounds[i]; + } + } + + //Momentum + offspring.Momentum[i] = Random.value * momentum + (offspring.Genes[i] - Gene); + } + + //Fitness + offspring.Fitness = Model.ComputeLoss(offspring.Genes); + } + + //Generates a random individual + private void Reroll(Individual individual) { + for(int i=0; i pool) { + double rankSum = (double)(PoolCount*(PoolCount+1)) / 2.0; + for(int i=0; i max) { + if(gene > max) { + gene = max + max - gene; + } + if(gene < min) { + gene = min + min - gene; + } + } + } + if(type == JointType.Continuous) { + //Overflow + while(gene < -PI || gene > PI) { + if(gene > PI) { + gene -= PI + PI; + } + if(gene < -PI) { + gene += PI + PI; + } + } + } + if(type == JointType.Floating) { + //nothing + } + return gene; + } + */ + + //Sorts the population by their fitness values (descending) + private void SortByFitness() { + System.Array.Sort(Population, + delegate(Individual a, Individual b) { + return a.Fitness.CompareTo(b.Fitness); + } + ); + } + + //In-place swap by references + private void Swap(ref Individual[] a, ref Individual[] b) { + Individual[] tmp = a; + a = b; + b = tmp; + } + + public Model GetModel() { + return Model; + } + + public double[] GetSolution() { + return Solution; + } + + public double GetFitness() { + return Fitness; + } + + public double[] GetLowerBounds() { + return LowerBounds; + } + + public double[] GetUpperBounds() { + return UpperBounds; + } + + public double[,] GetGeneLandscape() { + double[,] values = new double[PopulationSize, Dimensionality]; + for(int i=0; i Function; + public System.Func Gradient; + public int Dimensionality; + public double[] Solution; + public double Value; + public double[] LowerBounds; + public double[] UpperBounds; + + private const double stpmin = 1e-20; + private const double stpmax = 1e20; + private int Corrections = 1; + private double Factor = 1e+5; + private double Tolerance = 0.0; + private int IPrint = 101; + private int TotalSize = 0; + + private double F; + private double[] G; + + private bool[] LSave; + private int[] ISave; + private double[] DSave; + + private int[] IWA; + private int[] NBD; + private double[] Work; + + private Task _Task; + private Task _CSave; + + private double NewF; + private double[] NewG; + + public BFGS(int dimensionality, System.Func function, System.Func gradient) { + Dimensionality = dimensionality; + Function = function; + Gradient = gradient; + UpperBounds = new double[Dimensionality]; + LowerBounds = new double[Dimensionality]; + Solution = new double[Dimensionality]; + + TotalSize = 2 * Corrections * Dimensionality + 11 * Corrections * Corrections + 5 * Dimensionality + 8 * Corrections; + NBD = new int[Dimensionality]; + for(int i=0; i Function { get; set; } + public System.Func Gradient { get; set; } + public double[] Solution; + public double Value; + public double[] UpperBounds; + public double[] LowerBounds; + + private const double stpmin = 1e-20; + private const double stpmax = 1e20; + private int corrections = 1; + private double[] work; + private double factr = 1e+5; + private double pgtol = 0.0; + private int n; + private int m; + private Task task = Task.None; + private Task csave = Task.None; + private int totalSize; + private int iprint; + private bool[] lsave; + private int[] nbd; + private int[] iwa; + private int[] isave; + private double f; + private double[] l; + private double[] u; + private double[] g; + private double[] dsave; + private double newF; + private double[] newG; + + public BFGS(int numberOfVariables, System.Func function, System.Func gradient) { + NumberOfVariables = numberOfVariables; + Function = function; + Gradient = gradient; + Solution = new double[numberOfVariables]; + UpperBounds = new double[numberOfVariables]; + LowerBounds = new double[numberOfVariables]; + + n = NumberOfVariables; + m = corrections; + totalSize = 2 * m * n + 11 * m * m + 5 * n + 8 * m; + iprint = 101; + lsave = new bool[4]; + nbd = new int[n]; + iwa = new int[3*n]; + isave = new int[44]; + l = new double[n]; + u = new double[n]; + g = new double[n]; + dsave = new double[29]; + work = new double[totalSize]; + } + + private void Prepare(double[] values) { + for(int i=0; i 0)) + { + if (((nbd[(i - (1)) + _nbd_offset] <= 2) + && (x[(i - (1)) + _x_offset] <= l[(i - (1)) + _l_offset]))) + { + if ((x[(i - (1)) + _x_offset] < l[(i - (1)) + _l_offset])) + { + prjctd = true; + x[(i - (1)) + _x_offset] = l[(i - (1)) + _l_offset]; + } + nbdd = (nbdd + 1); + } + else if (((nbd[(i - (1)) + _nbd_offset] >= 2) + && (x[(i - (1)) + _x_offset] >= u[(i - (1)) + _u_offset]))) + { + if ((x[(i - (1)) + _x_offset] > u[(i - (1)) + _u_offset])) + { + prjctd = true; + x[(i - (1)) + _x_offset] = u[(i - (1)) + _u_offset]; + } + nbdd = (nbdd + 1); + } + } + } + } + // + // c Initialize iwhere and assign values to cnstnd and boxed. + // + { + for (i = 1; i <= n; i++) + { + if ((nbd[(i - (1)) + _nbd_offset] != 2)) + { + boxed = false; + } + if ((nbd[(i - (1)) + _nbd_offset] == 0)) + { + // this variable is always free + iwhere[(i - (1)) + _iwhere_offset] = -1; + // + // otherwise set x(i)=mid(x(i), u(i), l(i)). + } + else + { + cnstnd = true; + if (((nbd[(i - (1)) + _nbd_offset] == 2) && + ((u[(i - (1)) + _u_offset] - l[(i - (1)) + _l_offset]) <= 0.0))) + { + // this variable is always fixed + iwhere[(i - (1)) + _iwhere_offset] = 3; + } + else + { + iwhere[(i - (1)) + _iwhere_offset] = 0; + } + } + } + } + + // + if ((iprint >= 0)) + { + if (prjctd) + { + // DISPLAY: The initial X is infeasible. Restart with its projection."; + } + if ((!cnstnd)) + { + // DISPLAY: "This problem is unconstrained." + } + } + + if ((iprint > 0)) + { + // DISPLAY: 'At X0 ',i9,' variables are exactly at the bounds'" (nbdd) + } + } + + + // + // c====================== The end of dpofa =============================== + // + // c + // c + // c dtrsl solves systems of the form + // c + // c t * x = b + // c or + // c trans(t) * x = b + // c + // c where t is a triangular matrix of order n. here trans(t) + // c denotes the transpose of the matrix t. + // c + // c on entry + // c + // c t double precision(ldt,n) + // c t contains the matrix of the system. the zero + // c elements of the matrix are not referenced, and + // c the corresponding elements of the array can be + // c used to store other information. + // c + // c ldt integer + // c ldt is the leading dimension of the array t. + // c + // c n integer + // c n is the order of the system. + // c + // c b double precision(n). + // c b contains the right hand side of the system. + // c + // c job integer + // c job specifies what kind of system is to be solved. + // c if job is + // c + // c 00 solve t*x=b, t lower triangular, + // c 01 solve t*x=b, t upper triangular, + // c 10 solve trans(t)*x=b, t lower triangular, + // c 11 solve trans(t)*x=b, t upper triangular. + // c + // c on return + // c + // c b b contains the solution, if info .eq. 0. + // c otherwise b is unaltered. + // c + // c info integer + // c info contains zero if the system is nonsingular. + // c otherwise info contains the index of + // c the first zero diagonal element of t. + // c + // c linpack. this version dated 08/14/78 . + // c g. w. stewart, university of maryland, argonne national lab. + // c + // c subroutines and functions + // c + // c blas daxpy,ddot + // c fortran mod + // c + // c internal variables + // c + // c + // c begin block permitting ...exits to 150 + // c + // c check for zero diagonal elements. + // c + + private static void dtrsl(double[] t, int _t_offset, int ldt, int n, + double[] b, int _b_offset, int job, ref int info) + { + double temp = 0.0d; + int Case = 0; + int j = 0; + int jj = 0; + + { + for (info = 1; info <= n; info++) + { + // ......exit + if ((t[(info - (1)) + (info - (1)) * (ldt) + _t_offset] == 0.0e0)) + { + return; + } + } + } + + info = 0; + + // + // determine the task and go to it. + // + Case = 1; + + if (((job) % (10) != 0)) + { + Case = 2; + } + + if ((((job) % (100) / 10) != 0)) + { + Case = (Case + 2); + } + + { + int _cg_tmp = Case; + if (_cg_tmp == 1) + goto L20; + else if (_cg_tmp == 2) + goto L50; + else if (_cg_tmp == 3) + goto L80; + else if (_cg_tmp == 4) + goto L110; + } + + // + // solve t*x=b for t lower triangular + // + L20: + + b[(1 - (1)) + _b_offset] = (b[(1 - (1)) + _b_offset] + / t[(1 - (1)) + (1 - (1)) * (ldt) + _t_offset]); + + if ((n < 2)) + { + return; + } + + { + for (j = 2; j <= n; j++) + { + temp = (-(b[((j - 1) - (1)) + _b_offset])); + daxpy(((n - j) + 1), temp, t, (j - (1)) + ((j - 1) + - (1)) * (ldt) + _t_offset, 1, b, (j - (1)) + _b_offset, 1); + b[(j - (1)) + _b_offset] = (b[(j - (1)) + _b_offset] / t[(j - (1)) + + (j - (1)) * (ldt) + _t_offset]); + } + } + return; + + // + // solve t*x=b for t upper triangular. + // + L50: + + b[(n - (1)) + _b_offset] = (b[(n - (1)) + _b_offset] + / t[(n - (1)) + (n - (1)) * (ldt) + _t_offset]); + + if ((n < 2)) + { + return; + } + + { + for (jj = 2; jj <= n; jj++) + { + j = ((n - jj) + 1); + temp = (-(b[((j + 1) - (1)) + _b_offset])); + daxpy(j, temp, + t, (1 - (1)) + ((j + 1) - (1)) * (ldt) + _t_offset, 1, + b, (1 - (1)) + _b_offset, 1); + + b[(j - (1)) + _b_offset] = (b[(j - (1)) + _b_offset] + / t[(j - (1)) + (j - (1)) * (ldt) + _t_offset]); + } + } + + return; + // + // solve trans(t)*x=b for t lower triangular. + // + L80: + b[(n - (1)) + _b_offset] = (b[(n - (1)) + _b_offset] + / t[(n - (1)) + (n - (1)) * (ldt) + _t_offset]); + + if ((n < 2)) + { + return; + } + + { + for (jj = 2; jj <= n; jj++) + { + j = ((n - jj) + 1); + b[(j - (1)) + _b_offset] = (b[(j - (1)) + _b_offset] + - ddot((jj - 1), t, ((j + 1) - (1)) + (j - (1)) * (ldt) + _t_offset, + 1, b, ((j + 1) - (1)) + _b_offset, 1)); + + b[(j - (1)) + _b_offset] = (b[(j - (1)) + _b_offset] / t[(j - (1)) + + (j - (1)) * (ldt) + _t_offset]); + } + } + return; + + // + // solve trans(t)*x=b for t upper triangular. + // + L110: + b[(1 - (1)) + _b_offset] = (b[(1 - (1)) + _b_offset] + / t[(1 - (1)) + (1 - (1)) * (ldt) + _t_offset]); + + if ((n < 2)) + { + return; + } + + { + for (j = 2; j <= n; j++) + { + b[(j - (1)) + _b_offset] = (b[(j - (1)) + _b_offset] + - ddot((j - 1), t, (1 - (1)) + (j - (1)) * (ldt) + _t_offset, + 1, b, (1 - (1)) + _b_offset, 1)); + + b[(j - (1)) + _b_offset] = (b[(j - (1)) + _b_offset] + / t[(j - (1)) + (j - (1)) * (ldt) + _t_offset]); + } + } + } + + // + // c======================= The end of active ============================= + // + // + // + // c ************ + // c + // c Subroutine bmv + // c + // c This subroutine computes the product of the 2m x 2m middle matrix + // c in the compact L-BFGS formula of B and a 2m vector v; + // c it returns the product in p. + // c + // c m is an integer variable. + // c On entry m is the maximum number of variable metric corrections + + // c used to define the limited memory matrix. + // c On exit m is unchanged. + // c + // c sy is a double precision array of dimension m x m. + // c On entry sy specifies the matrix S'Y. + // c On exit sy is unchanged. + // c + // c wt is a double precision array of dimension m x m. + // c On entry wt specifies the upper triangular matrix J' which is + // c the Cholesky factor of (thetaS'S+LD^(-1)L'). + // c On exit wt is unchanged. + // c + // c col is an integer variable. + // c On entry col specifies the number of s-vectors (or y-vectors) + // c stored in the compact L-BFGS formula. + // c On exit col is unchanged. + // c + // c v is a double precision array of dimension 2col. + // c On entry v specifies vector v. + // c On exit v is unchanged. + // c + // c p is a double precision array of dimension 2col. + // c On entry p is unspecified. + // c On exit p is the product Mv. + // c + // c info is an integer variable. + // c On entry info is unspecified. + // c On exit info = 0 for normal return, + // c = nonzero for abnormal return when the system + // c to be solved by dtrsl is singular. + // c + // c Subprograms called: + // c + // c Linpack ... + // c + // c + // c * * * + // c + // c NEOS, November 1994. (Latest revision June 1996.) + // c Optimization Technology Center. + // c Argonne National Laboratory and Northwestern University. + // c Written by + // c Ciyou Zhu + // c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. + // c + // c + // c ************ + // + // + private static void bmv(int m, double[] sy, int _sy_offset, double[] wt, int _wt_offset, + int col, double[] v, int _v_offset, double[] p, int _p_offset, ref int info) + { + + int i = 0; + int k = 0; + int i2 = 0; + double sum = 0.0d; + + if ((col == 0)) + return; + + // + // PART I: solve [ D^(1/2) O ] [ p1 ] = [ v1 ] + // [ -L*D^(-1/2) J ] [ p2 ] [ v2 ]. + // + // solve Jp2=v2+LD^(-1)v1. + // + p[((col + 1) - (1)) + _p_offset] = v[((col + 1) - (1)) + _v_offset]; + { + for (i = 2; i <= col; i++) + { + i2 = (col + i); + sum = 0.0e0; + { + for (k = 1; k <= (i - 1); k++) + { + sum = (sum + ((sy[(i - (1)) + (k - (1)) * (m) + + _sy_offset] * v[(k - (1)) + _v_offset]) / sy[(k - (1)) + + (k - (1)) * (m) + _sy_offset])); + } + } + + p[(i2 - (1)) + _p_offset] = (v[(i2 - (1)) + _v_offset] + sum); + } + } + + // Solve the triangular system + dtrsl(wt, _wt_offset, m, col, p, + ((col + 1) - (1)) + _p_offset, 11, ref info); + + // + // solve D^(1/2)p1=v1. + { + for (i = 1; i <= col; i++) + { + p[(i - (1)) + _p_offset] = (v[(i - (1)) + _v_offset] + / System.Math.Sqrt(sy[(i - (1)) + (i - (1)) * (m) + _sy_offset])); + } + } + + // + // PART II: solve [ -D^(1/2) D^(-1/2)*L' ] [ p1 ] = [ p1 ] + // [ 0 J' ] [ p2 ] [ p2 ]. + // + // solve J^Tp2=p2. + // + dtrsl(wt, _wt_offset, m, col, + p, ((col + 1) - (1)) + _p_offset, 01, ref info); + + // + // compute p1 = -D^(-1/2)(p1-D^(-1/2)L'p2) + // = -D^(-1/2)p1+D^(-1)L'p2. + { + for (i = 1; i <= col; i++) + { + p[(i - (1)) + _p_offset] = (-((p[(i - (1)) + _p_offset] + / System.Math.Sqrt(sy[(i - (1)) + (i - (1)) * (m) + _sy_offset])))); + } + } + + { + for (i = 1; i <= col; i++) + { + sum = 0.0e0; + { + for (k = (i + 1); k <= col; k++) + { + sum = (sum + ((sy[(k - (1)) + (i - (1)) * (m) + _sy_offset] + * p[((col + k) - (1)) + _p_offset]) / sy[(i - (1)) + (i - (1)) * (m) + _sy_offset])); + } + } + p[(i - (1)) + _p_offset] = (p[(i - (1)) + _p_offset] + sum); + } + } + } + + // + // c======================== The end of bmv =============================== + // + // + // c ************ + // c + // c Subroutine cauchy + // c + // c For given x, l, u, g (with sbgnrm > 0), and a limited memory + // c BFGS matrix B defined in terms of matrices WY, WS, WT, and + // c scalars head, col, and theta, this subroutine computes the + // c generalized Cauchy point (GCP), defined as the first local + // c minimizer of the quadratic + // c + // c Q(x + s) = g's + 1/2 s'Bs + // c + // c along the projected gradient direction P(x-tg,l,u). + // c The routine returns the GCP in xcp. + // c + // c n is an integer variable. + // c On entry n is the dimension of the problem. + // c On exit n is unchanged. + // c + // c x is a double precision array of dimension n. + // c On entry x is the starting point for the GCP computation. + // c On exit x is unchanged. + // c + // c l is a double precision array of dimension n. + // c On entry l is the lower bound of x. + // c On exit l is unchanged. + // c + // c u is a double precision array of dimension n. + // c On entry u is the upper bound of x. + // c On exit u is unchanged. + // c + // c nbd is an integer array of dimension n. + // c On entry nbd represents the type of bounds imposed on the + // c variables, and must be specified as follows: + // c nbd(i)=0 if x(i) is unbounded, + // c 1 if x(i) has only a lower bound, + // c 2 if x(i) has both lower and upper bounds, and + // c 3 if x(i) has only an upper bound. + // c On exit nbd is unchanged. + // c + // c g is a double precision array of dimension n. + // c On entry g is the gradient of f(x). g must be a nonzero vector. + // c On exit g is unchanged. + // c + // c iorder is an integer working array of dimension n. + // c iorder will be used to store the breakpoints in the piecewise + // c linear path and free variables encountered. On exit, + // c iorder(1),...,iorder(nleft) are indices of breakpoints + // c which have not been encountered; + // c iorder(nleft+1),...,iorder(nbreak) are indices of + // c encountered breakpoints; and + // c iorder(nfree),...,iorder(n) are indices of variables which + // c have no bound constraits along the search direction. + // c + // c iwhere is an integer array of dimension n. + // c On entry iwhere indicates only the permanently fixed (iwhere=3) + + // c or free (iwhere= -1) components of x. + // c On exit iwhere records the status of the current x variables. + // c iwhere(i)=-3 if x(i) is free and has bounds, but is not moved + // c 0 if x(i) is free and has bounds, and is moved + // c 1 if x(i) is fixed at l(i), and l(i) .ne. u(i) + // c 2 if x(i) is fixed at u(i), and u(i) .ne. l(i) + // c 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i) + // c -1 if x(i) is always free, i.e., it has no bounds. + // c + // c t is a double precision working array of dimension n. + // c t will be used to store the break points. + // c + // c d is a double precision array of dimension n used to store + // c the Cauchy direction P(x-tg)-x. + // c + // c xcp is a double precision array of dimension n used to return the + + // c GCP on exit. + // c + // c m is an integer variable. + // c On entry m is the maximum number of variable metric corrections + // c used to define the limited memory matrix. + // c On exit m is unchanged. + // c + // c ws, wy, sy, and wt are double precision arrays. + // c On entry they store information that defines the + // c limited memory BFGS matrix: + // c ws(n,m) stores S, a set of s-vectors; + // c wy(n,m) stores Y, a set of y-vectors; + // c sy(m,m) stores S'Y; + // c wt(m,m) stores the + // c Cholesky factorization of (theta*S'S+LD^(-1)L'). + // c On exit these arrays are unchanged. + // c + // c theta is a double precision variable. + // c On entry theta is the scaling factor specifying B_0 = theta I. + // c On exit theta is unchanged. + // c + // c col is an integer variable. + // c On entry col is the actual number of variable metric + // c corrections stored so far. + // c On exit col is unchanged. + // c + // c head is an integer variable. + // c On entry head is the location of the first s-vector (or y-vector + // c in S (or Y). + // c On exit col is unchanged. + // c + // c p is a double precision working array of dimension 2m. + // c p will be used to store the vector p = W^(T)d. + // c + // c c is a double precision working array of dimension 2m. + // c c will be used to store the vector c = W^(T)(xcp-x). + // c + // c wbp is a double precision working array of dimension 2m. + // c wbp will be used to store the row of W corresponding + // c to a breakpoint. + // c + // c v is a double precision working array of dimension 2m. + // c + // c nseg is an integer variable. + // c On exit nseg records the number of quadratic segments explored + // c in searching for the GCP. + // c + // c sg and yg are double precision arrays of dimension m. + // c On entry sg and yg store S'g and Y'g correspondingly. + // c On exit they are unchanged. + // c + // c iprint is an INTEGER variable that must be set by the user. + // c It controls the frequency and type of output generated: + // c iprint<0 no output is generated; + // c iprint=0 print only one line at the last iteration; + // c 0100 print details of every iteration including x and g; + // c When iprint > 0, the file iterate.dat will be created to + // c summarize the iteration. + // c + // c sbgnrm is a double precision variable. + // c On entry sbgnrm is the norm of the projected gradient at x. + // c On exit sbgnrm is unchanged. + // c + // c info is an integer variable. + // c On entry info is 0. + // c On exit info = 0 for normal return, + // c = nonzero for abnormal return when the the system + // c used in routine bmv is singular. + // c + // c Subprograms called: + // c + // c L-BFGS-B Library ... hpsolb, bmv. + // c + // c Linpack ... dscal dcopy, + // c + // c + // c References: + // c + // c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited + // c memory algorithm for bound constrained optimization'', + // c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. + // c + // c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN + // c Subroutines for Large Scale Bound Constrained Optimization'' + // c Tech. Report, NAM-11, EECS Department, Northwestern University, + + // c 1994. + // c + // c (Postscript files of these papers are available via anonymous + // c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) + // c + // c * * * + // c + // c NEOS, November 1994. (Latest revision June 1996.) + // c Optimization Technology Center. + // c Argonne National Laboratory and Northwestern University. + // c Written by + // c Ciyou Zhu + // c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. + // c + // c + // c ************ + // + // + // c Check the status of the variables, reset iwhere(i) if necessary; + // c compute the Cauchy direction d and the breakpoints t; initialize + // c the derivative f1 and the vector p = W'd (for theta = 1). + // + private static void cauchy(int n, double[] x, int _x_offset, double[] l, int _l_offset, + double[] u, int _u_offset, int[] nbd, int _nbd_offset, double[] g, int _g_offset, + int[] iorder, int _iorder_offset, int[] iwhere, int _iwhere_offset, double[] t, int _t_offset, + double[] d, int _d_offset, double[] xcp, int _xcp_offset, int m, + double[] wy, int _wy_offset, double[] ws, int _ws_offset, double[] sy, int _sy_offset, + double[] wt, int _wt_offset, double theta, int col, + int head, double[] p, int _p_offset, double[] c, int _c_offset, + double[] wbp, int _wbp_offset, double[] v, int _v_offset, ref int nseg, + int iprint, double sbgnrm, ref int info, double epsmch) + { + + bool xlower = false; + bool xupper = false; + bool bnded = false; + int i = 0; + int j = 0; + int col2 = 0; + int nfree = 0; + int nbreak = 0; + int pointr = 0; + int ibp = 0; + int nleft = 0; + int ibkmin = 0; + int iter = 0; + double f1 = 0.0d; + double f2 = 0.0d; + double dt = 0.0d; + double dtm = 0.0d; + double tsum = 0.0d; + double dibp = 0.0d; + double zibp = 0.0d; + double dibp2 = 0.0d; + double bkmin = 0.0d; + double tu = 0.0d; + double tl = 0.0d; + double wmc = 0.0d; + double wmp = 0.0d; + double wmw = 0.0d; + double tj = 0.0d; + double tj0 = 0.0d; + double neggi = 0.0d; + double f2_org = 0.0d; + + if ((sbgnrm <= 0.0)) + { + if ((iprint >= 0)) + { + // DISPLAY: Subgnorm = 0. GCP = X. + } + + dcopy(n, x, _x_offset, 1, xcp, _xcp_offset, 1); + return; + } + + bnded = true; + nfree = (n + 1); + nbreak = 0; + ibkmin = 0; + bkmin = 0.0; + col2 = (2 * col); + f1 = 0.0; + if ((iprint >= 99)) + { + // DISPLAY: '---------------- CAUCHY entered-------------------' + } + // + // c We set p to zero and build it up as we determine d. + // + { + for (i = 1; i <= col2; i++) + { + p[(i - (1)) + _p_offset] = 0.0; + } + } + // + // c In the following loop we determine for each variable its bound + // c status and its breakpoint, and update p accordingly. + // c Smallest breakpoint is identified. + // + { + for (i = 1; i <= n; i++) + { + neggi = (-(g[(i - (1)) + _g_offset])); + + if (((iwhere[(i - (1)) + _iwhere_offset] != 3) + && (iwhere[(i - (1)) + _iwhere_offset] != -1))) + { + // c if x(i) is not a constant and has bounds, + // c compute the difference between x(i) and its bounds. + if ((nbd[(i - (1)) + _nbd_offset] <= 2)) + { + tl = (x[(i - (1)) + _x_offset] - l[(i - (1)) + _l_offset]); + } + if ((nbd[(i - (1)) + _nbd_offset] >= 2)) + { + tu = (u[(i - (1)) + _u_offset] - x[(i - (1)) + _x_offset]); + } + // + // c If a variable is close enough to a bound + // c we treat it as at bound. + xlower = ((nbd[(i - (1)) + _nbd_offset] <= 2) && (tl <= 0.0)); + xupper = ((nbd[(i - (1)) + _nbd_offset] >= 2) && (tu <= 0.0)); + // + // c reset iwhere(i). + iwhere[(i - (1)) + _iwhere_offset] = 0; + if (xlower) + { + if ((neggi <= 0.0)) + { + iwhere[(i - (1)) + _iwhere_offset] = 1; + } + } + else if (xupper) + { + if ((neggi >= 0.0)) + { + iwhere[(i - (1)) + _iwhere_offset] = 2; + } + } + else + { + if ((System.Math.Abs(neggi) <= 0.0)) + { + iwhere[(i - (1)) + _iwhere_offset] = -3; + } + } + } + pointr = head; + if (((iwhere[(i - (1)) + _iwhere_offset] != 0) && (iwhere[(i - (1)) + _iwhere_offset] != -1))) + { + d[(i - (1)) + _d_offset] = 0.0; + } + else + { + d[(i - (1)) + _d_offset] = neggi; + f1 = (f1 - (neggi * neggi)); + + // calculate p := p - W'e_i* (g_i). + { + for (j = 1; j <= col; j++) + { + p[(j - (1)) + _p_offset] = (p[(j - (1)) + + _p_offset] + (wy[(i - (1)) + (pointr - (1)) * (n) + _wy_offset] * neggi)); + p[((col + j) - (1)) + _p_offset] = (p[((col + j) - (1)) + + _p_offset] + (ws[(i - (1)) + (pointr - (1)) * (n) + _ws_offset] * neggi)); + pointr = ((pointr) % (m) + 1); + } + } + if ((((nbd[(i - (1)) + _nbd_offset] <= 2) + && (nbd[(i - (1)) + _nbd_offset] != 0)) && (neggi < 0.0))) + { + // x(i) + d(i) is bounded; compute t(i). + + nbreak = (nbreak + 1); + iorder[(nbreak - (1)) + _iorder_offset] = i; + t[(nbreak - (1)) + _t_offset] = (tl / ((-(neggi)))); + if (((nbreak == 1) || (t[(nbreak - (1)) + _t_offset] < bkmin))) + { + bkmin = t[(nbreak - (1)) + _t_offset]; + ibkmin = nbreak; + } + } + else if (((nbd[(i - (1)) + _nbd_offset] >= 2) && (neggi > 0.0))) + { + // x(i) + d(i) is bounded; compute t(i). + + nbreak = (nbreak + 1); + iorder[(nbreak - (1)) + _iorder_offset] = i; + t[(nbreak - (1)) + _t_offset] = (tu / neggi); + if (((nbreak == 1) || (t[(nbreak - (1)) + _t_offset] < bkmin))) + { + bkmin = t[(nbreak - (1)) + _t_offset]; + ibkmin = nbreak; + } + } + else + { + // x(i) + d(i) is not bounded. + nfree = (nfree - 1); + iorder[(nfree - (1)) + _iorder_offset] = i; + if ((System.Math.Abs(neggi) > 0.0)) + { + bnded = false; + } + } + } + } + } + // + // The indices of the nonzero components of d are now stored + // in iorder(1),...,iorder(nbreak) and iorder(nfree),...,iorder(n). + // The smallest of the nbreak breakpoints is in t(ibkmin)=bkmin. + // + if ((theta != 1.0)) + { + // complete the initialization of p for theta not= one. + dscal(col, theta, p, ((col + 1) - (1)) + _p_offset, 1); + } + + // + // c Initialize GCP xcp = x. + // + dcopy(n, x, _x_offset, 1, xcp, _xcp_offset, 1); + + if (((nbreak == 0) && (nfree == (n + 1)))) + { + // is a zero vector, return with the initial xcp as GCP. + return; + } + + // + // c Initialize c = W'(xcp - x) = 0. + // + { + for (j = 1; j <= col2; j++) + { + c[(j - (1)) + _c_offset] = 0.0; + } + } + + // + // c Initialize derivative f2. + // + f2 = (-((theta * f1))); + f2_org = f2; + if ((col > 0)) + { + bmv(m, sy, _sy_offset, wt, _wt_offset, + col, p, _p_offset, v, _v_offset, ref info); + + f2 = (f2 - BFGS.ddot(col2, v, _v_offset, 1, p, _p_offset, 1)); + } + + dtm = (-((f1 / f2))); + tsum = 0.0; + nseg = 1; + + if ((iprint >= 99)) + { + // DISPLAY: "There are " + nbreak + " breakpoints " + } + + // + // c If there are no breakpoints, locate the GCP and return. + // + if ((nbreak == 0)) + goto L888; + + nleft = nbreak; + iter = 1; + tj = 0.0; + + // + // c------------------- the beginning of the loop ------------------------- + // + L777: + // + // Find the next smallest breakpoint; + // compute dt = t(nleft) - t(nleft + 1). + // + tj0 = tj; + + if ((iter == 1)) + { + // Since we already have the smallest breakpoint we need not do + // heapsort yet. Often only one breakpoint is used and the + // cost of heapsort is avoided. + tj = bkmin; + ibp = iorder[(ibkmin - (1)) + _iorder_offset]; + } + else + { + if ((iter == 2)) + { + // Replace the already used smallest breakpoint with the + // breakpoint numbered nbreak > nlast, before heapsort call. + + if ((ibkmin != nbreak)) + { + t[(ibkmin - (1)) + _t_offset] = t[(nbreak - (1)) + _t_offset]; + iorder[(ibkmin - (1)) + _iorder_offset] = iorder[(nbreak - (1)) + _iorder_offset]; + } + // Update heap structure of breakpoints + // (if iter=2, initialize heap). + } + hpsolb(nleft, t, _t_offset, iorder, _iorder_offset, (iter - 2)); + tj = t[(nleft - (1)) + _t_offset]; + ibp = iorder[(nleft - (1)) + _iorder_offset]; + } + // + dt = (tj - tj0); + // + if (((dt != 0.0) && (iprint >= 100))) + { + // DISPLAY: nseg, f1, f2 + // "/,'Piece ',i3,' --f1, f2 at start point ',1p,2(1x,d11.4)" + // + // DISPLAY: dt, + // "'Distance to the next break point = ',1p,d11.4" + // + // DISPLAY: dtm, + // "'Distance to the stationary point = ',1p,d11.4" + } + + // + // If a minimizer is within this interval, locate the GCP and return. + // + if ((dtm < dt)) + goto L888; + + // + // Otherwise fix one variable and + // reset the corresponding component of d to zero. + // + tsum = (tsum + dt); + nleft = (nleft - 1); + iter = (iter + 1); + dibp = d[(ibp - (1)) + _d_offset]; + d[(ibp - (1)) + _d_offset] = 0.0; + + if ((dibp > 0.0)) + { + zibp = (u[(ibp - (1)) + _u_offset] - x[(ibp - (1)) + _x_offset]); + xcp[(ibp - (1)) + _xcp_offset] = u[(ibp - (1)) + _u_offset]; + iwhere[(ibp - (1)) + _iwhere_offset] = 2; + } + else + { + zibp = (l[(ibp - (1)) + _l_offset] - x[(ibp - (1)) + _x_offset]); + xcp[(ibp - (1)) + _xcp_offset] = l[(ibp - (1)) + _l_offset]; + iwhere[(ibp - (1)) + _iwhere_offset] = 1; + } + + if ((iprint >= 100)) + { + // DISPLAY: "Variable " + ibp + " is fixed." + } + + if (((nleft == 0) && (nbreak == n))) + { + // all n variables are fixed, + // return with xcp as GCP. + dtm = dt; + goto L999; + } + + // + // Update the derivative information. + // + nseg = (nseg + 1); + dibp2 = (System.Math.Pow(dibp, 2)); + + // + // Update f1 and f2. + // + // temporarily set f1 and f2 for col=0. + // + f1 = (((f1 + (dt * f2)) + dibp2) - ((theta * dibp) * zibp)); + f2 = (f2 - (theta * dibp2)); + // + if ((col > 0)) + { + // update c = c + dt*p. + daxpy(col2, dt, p, _p_offset, 1, c, _c_offset, 1); + + // choose wbp, + // the row of W corresponding to the breakpoint encountered. + pointr = head; + { + for (j = 1; j <= col; j++) + { + wbp[(j - (1)) + _wbp_offset] = wy[(ibp - (1)) + + (pointr - (1)) * (n) + _wy_offset]; + + wbp[((col + j) - (1)) + _wbp_offset] = (theta * ws[(ibp + - (1)) + (pointr - (1)) * (n) + _ws_offset]); + + pointr = ((pointr) % (m) + 1); + } + } + + // compute (wbp)Mc, (wbp)Mp, and (wbp)M(wbp)'. + bmv(m, sy, _sy_offset, wt, _wt_offset, col, wbp, + _wbp_offset, v, _v_offset, ref info); + + wmc = ddot(col2, c, _c_offset, 1, v, _v_offset, 1); + wmp = ddot(col2, p, _p_offset, 1, v, _v_offset, 1); + wmw = ddot(col2, wbp, _wbp_offset, 1, v, _v_offset, 1); + + // update p = p - dibp*wbp. + daxpy(col2, (-(dibp)), wbp, _wbp_offset, 1, p, _p_offset, 1); + + // complete updating f1 and f2 while col > 0. + f1 = (f1 + (dibp * wmc)); + f2 = ((f2 + ((2.0e0 * dibp) * wmp)) - (dibp2 * wmw)); + } + + + f2 = System.Math.Max((epsmch * f2_org), f2); + + if ((nleft > 0)) + { + dtm = (-((f1 / f2))); + goto L777; + // to repeat the loop for unsearched intervals. + } + else if (bnded) + { + f1 = 0.0; + f2 = 0.0; + dtm = 0.0; + } + else + { + dtm = (-((f1 / f2))); + } + + // + // c------------------- the end of the loop ------------------------------- + // + L888: + + if ((iprint >= 99)) + { + // DISPLAY: "GCP found in this segment", + // + // DISPLAY: nseg, f1, f2 + // 'Piece ',i3,' --f1, f2 at start point ',1p,2(1x,d11.4)" + // dtm + // "'Distance to the stationary point = ',1p,d11.4" + } + if ((dtm <= 0.0)) + { + dtm = 0.0; + } + tsum = (tsum + dtm); + + // + // Move free variables (i.e., the ones w/o breakpoints) and + // the variables whose breakpoints haven't been reached. + // + daxpy(n, tsum, d, _d_offset, 1, xcp, _xcp_offset, 1); + // + L999: + // + // Update c = c + dtm*p = W'(x^c - x) + // which will be used in computing r = Z'(B(x^c - x) + g). + // + if ((col > 0)) + { + daxpy(col2, dtm, p, _p_offset, 1, c, _c_offset, 1); + } + + + if ((iprint > 100)) + { + for (i = 1; i <= n; i++) + { + // DISPLAY: xcp[(i - (1)) + _xcp_offset]; + } + // DISPLAY: "'Cauchy X = ',/,(4x,1p,6(1x,d11.4))" + } + + if ((iprint >= 99)) + { + // DISPLAY: '---------------- exit CAUCHY----------------------' + } + + } + + // + // c====================== The end of cauchy ============================== + // + // + // + // c ************ + // c + // c Subroutine cmprlb + // c + // c This subroutine computes r=-Z'B(xcp-xk)-Z'g by using + // c wa(2m+1)=W'(xcp-x) from subroutine cauchy. + // c + // c Subprograms called: + // c + // c L-BFGS-B Library ... bmv. + // c + // c + // c * * * + // c + // c NEOS, November 1994. (Latest revision June 1996.) + // c Optimization Technology Center. + // c Argonne National Laboratory and Northwestern University. + // c Written by + // c Ciyou Zhu + // c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. + // c + // c + // c ************ + // + // + private static void cmprlb(int n, int m, double[] x, int _x_offset, double[] g, int _g_offset, double[] ws, int _ws_offset, double[] wy, int _wy_offset, double[] sy, int _sy_offset, + double[] wt, int _wt_offset, double[] z, int _z_offset, double[] r, int _r_offset, + double[] wa, int _wa_offset, int[] index, int _index_offset, double theta, + int col, int head, int nfree, bool cnstnd, ref int info) + { + + int i = 0; + int j = 0; + int k = 0; + int pointr = 0; + double a1 = 0.0d; + double a2 = 0.0d; + + if (((!cnstnd) && (col > 0))) + { + { + for (i = 1; i <= n; i++) + { + r[(i - (1)) + _r_offset] = (-(g[(i - (1)) + _g_offset])); + } + } + } + else + { + { + for (i = 1; i <= nfree; i++) + { + k = index[(i - (1)) + _index_offset]; + r[(i - (1)) + _r_offset] = ((-((theta * ((z[(k - (1)) + + _z_offset] - x[(k - (1)) + _x_offset]))))) - g[(k - (1)) + _g_offset]); + } + } + + bmv(m, sy, _sy_offset, wt, _wt_offset, col, wa, + (((2 * m) + 1) - (1)) + _wa_offset, wa, (1 - (1)) + _wa_offset, ref info); + + pointr = head; + + { + for (j = 1; j <= col; j++) + { + a1 = wa[(j - (1)) + _wa_offset]; + a2 = (theta * wa[((col + j) - (1)) + _wa_offset]); + { + for (i = 1; i <= nfree; i++) + { + k = index[(i - (1)) + _index_offset]; + r[(i - (1)) + _r_offset] = ((r[(i - (1)) + _r_offset] + + (wy[(k - (1)) + (pointr - (1)) * (n) + _wy_offset] * a1)) + + (ws[(k - (1)) + (pointr - (1)) * (n) + _ws_offset] * a2)); + } + } + + pointr = ((pointr) % (m) + 1); + } + } + } + } + + // c====================== The end of subsm =============================== + // + // c ********** + // c + // c Subroutine dcsrch + // c + // c This subroutine finds a step that satisfies a sufficient + // c decrease condition and a curvature condition. + // c + // c Each call of the subroutine updates an interval with + // c endpoints stx and sty. The interval is initially chosen + // c so that it contains a minimizer of the modified function + // c + // c psi(stp) = f(stp) - f(0) - ftol*stp*f'(0). + // c + // c If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the + // c interval is chosen so that it contains a minimizer of f. + // c + // c The algorithm is designed to find a step that satisfies + // c the sufficient decrease condition + // c + // c f(stp) <= f(0) + ftol*stp*f'(0), + // c + // c and the curvature condition + // c + // c abs(f'(stp)) <= gtol*abs(f'(0)). + // c + // c If ftol is less than gtol and if, for example, the function + // c is bounded below, then there is always a step which satisfies + // c both conditions. + // c + // c If no step can be found that satisfies both conditions, then + // c the algorithm stops with a warning. In this case stp only + // c satisfies the sufficient decrease condition. + // c + // c A typical invocation of dcsrch has the following outline: + // c + // c task = 'START' + // c 10 continue + // c call dcsrch( ... ) + // c if (task .eq. 'FG') then + // c Evaluate the function and the gradient at stp + // c goto 10 + // c end if + // c + // c NOTE: The user must no alter work arrays between calls. + // c + // c The subroutine statement is + // c + // c subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax, + // c task,isave,dsave) + // c where + // c + // c f is a double precision variable. + // c On initial entry f is the value of the function at 0. + // c On subsequent entries f is the value of the + // c function at stp. + // c On exit f is the value of the function at stp. + // c + // c g is a double precision variable. + // c On initial entry g is the derivative of the function at 0. + // c On subsequent entries g is the derivative of the + // c function at stp. + // c On exit g is the derivative of the function at stp. + // c + // c stp is a double precision variable. + // c On entry stp is the current estimate of a satisfactory + // c step. On initial entry, a positive initial estimate + // c must be provided. + // c On exit stp is the current estimate of a satisfactory step + // c if task = 'FG'. If task = 'CONV' then stp satisfies + // c the sufficient decrease and curvature condition. + // c + // c ftol is a double precision variable. + // c On entry ftol specifies a nonnegative tolerance for the + // c sufficient decrease condition. + // c On exit ftol is unchanged. + // c + // c gtol is a double precision variable. + // c On entry gtol specifies a nonnegative tolerance for the + // c curvature condition. + // c On exit gtol is unchanged. + // c + // c xtol is a double precision variable. + // c On entry xtol specifies a nonnegative relative tolerance + // c for an acceptable step. The subroutine exits with a + // c warning if the relative difference between sty and stx + // c is less than xtol. + // c On exit xtol is unchanged. + // c + // c stpmin is a double precision variable. + // c On entry stpmin is a nonnegative lower bound for the step. + // c On exit stpmin is unchanged. + // c + // c stpmax is a double precision variable. + // c On entry stpmax is a nonnegative upper bound for the step. + // c On exit stpmax is unchanged. + // c + // c task is a character variable of length at least 60. + // c On initial entry task must be set to 'START'. + // c On exit task indicates the required action: + // c + // c If task(1:2) = 'FG' then evaluate the function and + // c derivative at stp and call dcsrch again. + // c + // c If task(1:4) = 'CONV' then the search is successful. + // c + // c If task(1:4) = 'WARN' then the subroutine is not able + // c to satisfy the convergence conditions. The exit value of + // c stp contains the best point found during the search. + // c + // c If task(1:5) = 'ERROR' then there is an error in the + // c input arguments. + // c + // c On exit with convergence, a warning or an error, the + // c variable task contains additional information. + // c + // c isave is an integer work array of dimension 2. + // c + // c dsave is a double precision work array of dimension 13. + // c + // c Subprograms called + // c + // c MINPACK-2 ... dcstep + // c + // c MINPACK-1 Project. June 1983. + // c Argonne National Laboratory. + // c Jorge J. More' and David J. Thuente. + // c + // c MINPACK-2 Project. October 1993. + // c Argonne National Laboratory and University of Minnesota. + // c Brett M. Averick, Richard G. Carter, and Jorge J. More'. + // c + // c ********** + // + // + // c Initialization block. + // + private static void dcsrch(double f, double g, ref double stp, double ftol, + double gtol, double xtol, double stpmin, double stpmax, ref Task task, + int[] isave, int _isave_offset, double[] dsave, int _dsave_offset) + { + bool brackt = false; + int stage = 0; + double finit = 0.0d; + double ftest = 0.0d; + double fm = 0.0d; + double fx = 0.0d; + double fxm = 0.0d; + double fy = 0.0d; + double fym = 0.0d; + double ginit = 0.0d; + double gtest = 0.0d; + double gm = 0.0d; + double gx = 0.0d; + double gxm = 0.0d; + double gy = 0.0d; + double gym = 0.0d; + double stx = 0.0d; + double sty = 0.0d; + double stmin = 0.0d; + double stmax = 0.0d; + double width = 0.0d; + double width1 = 0.0d; + + if (task == Task.Start) + { + // + // Initialize local variables. + // + brackt = false; + stage = 1; + finit = f; + ginit = g; + gtest = (ftol * ginit); + width = (stpmax - stpmin); + width1 = (width / 0.5); + + // + // The variables stx, fx, gx contain the values of the step, + // function, and derivative at the best step. + // The variables sty, fy, gy contain the value of the step, + // function, and derivative at sty. + // The variables stp, f, g contain the values of the step, + // function, and derivative at stp. + // + stx = 0.0; + fx = finit; + gx = ginit; + sty = 0.0; + fy = finit; + gy = ginit; + stmin = 0.0; + stmax = (stp + (4.0 * stp)); + task = Task.FG; + goto L1000; + } + else + { + // Restore local variables. + if ((isave[(1 - (1)) + _isave_offset] == 1)) + { + brackt = true; + } + else + { + brackt = false; + } + stage = isave[(2 - (1)) + _isave_offset]; + ginit = dsave[(1 - (1)) + _dsave_offset]; + gtest = dsave[(2 - (1)) + _dsave_offset]; + gx = dsave[(3 - (1)) + _dsave_offset]; + gy = dsave[(4 - (1)) + _dsave_offset]; + finit = dsave[(5 - (1)) + _dsave_offset]; + fx = dsave[(6 - (1)) + _dsave_offset]; + fy = dsave[(7 - (1)) + _dsave_offset]; + stx = dsave[(8 - (1)) + _dsave_offset]; + sty = dsave[(9 - (1)) + _dsave_offset]; + stmin = dsave[(10 - (1)) + _dsave_offset]; + stmax = dsave[(11 - (1)) + _dsave_offset]; + width = dsave[(12 - (1)) + _dsave_offset]; + width1 = dsave[(13 - (1)) + _dsave_offset]; + } + + // + // c If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the + // c algorithm enters the second stage. + // + ftest = (finit + (stp * gtest)); + if ((((stage == 1) && (f <= ftest)) && (g >= 0.0))) + { + stage = 2; + } + + // + // c Test for warnings. + // + if ((brackt && (((stp <= stmin) || (stp >= stmax))))) + { + task = Task.Warning; + } + if ((brackt && ((stmax - stmin) <= (xtol * stmax)))) + { + task = Task.Warning; + } + if ((((stp == stpmax) && (f <= ftest)) && (g <= gtest))) + { + task = Task.Warning; + } + if (((stp == stpmin) && (((f > ftest) || (g >= gtest))))) + { + task = Task.Warning; + } + + // + // c Test for convergence. + // + if (((f <= ftest) && (System.Math.Abs(g) <= (gtol * ((-(ginit))))))) + { + task = Task.Convergence; + } + // + // c Test for termination. + // + if (task == Task.Warning || task == Task.Convergence) + { + goto L1000; + } + + // + // c A modified function is used to predict the step during the + // c first stage if a lower function value has been obtained but + // c the decrease is not sufficient. + // + if ((((stage == 1) && (f <= fx)) && (f > ftest))) + { + // + // c Define the modified function and derivative values. + // + fm = (f - (stp * gtest)); + fxm = (fx - (stx * gtest)); + fym = (fy - (sty * gtest)); + gm = (g - gtest); + gxm = (gx - gtest); + gym = (gy - gtest); + + // + // Call dcstep to update stx, sty, and to compute the new step. + // + dcstep(ref stx, ref fxm, ref gxm, ref sty, ref fym, ref gym, + ref stp, fm, gm, ref brackt, stmin, stmax); + + // + // Reset the function and derivative values for f. + // + fx = (fxm + (stx * gtest)); + fy = (fym + (sty * gtest)); + gx = (gxm + gtest); + gy = (gym + gtest); + // + } + else + { + // + // Call dcstep to update stx, sty, and to compute the new step. + // + dcstep(ref stx, ref fx, ref gx, ref sty, ref fy, ref gy, + ref stp, f, g, ref brackt, stmin, stmax); + } + + // + // c Decide if a bisection step is needed. + // + if (brackt) + { + if ((System.Math.Abs((sty - stx)) >= (0.6600000000000000310862446895043831318617 * width1))) + { + stp = (stx + (0.5 * ((sty - stx)))); + } + width1 = width; + width = System.Math.Abs((sty - stx)); + } + // + // c Set the minimum and maximum steps allowed for stp. + // + if (brackt) + { + stmin = System.Math.Min(stx, sty); + stmax = System.Math.Max(stx, sty); + } + else + { + stmin = (stp + (1.100000000000000088817841970012523233891 * ((stp - stx)))); + stmax = (stp + (4.0 * ((stp - stx)))); + } + + // + // c Force the step to be within the bounds stpmax and stpmin. + // + stp = System.Math.Max(stp, stpmin); + stp = System.Math.Min(stp, stpmax); + // + // c If further progress is not possible, let stp be the best + // c point obtained during the search. + // + if (((brackt && (((stp <= stmin) || (stp >= stmax)))) || ((brackt && ((stmax - stmin) <= (xtol * stmax)))))) + { + stp = stx; + } + // + // c Obtain another function and derivative. + // + task = Task.FG; + + // + L1000: + + // + // c Save local variables. + // + if (brackt) + { + isave[(1 - (1)) + _isave_offset] = 1; + } + else + { + isave[(1 - (1)) + _isave_offset] = 0; + } + + isave[(2 - (1)) + _isave_offset] = stage; + dsave[(1 - (1)) + _dsave_offset] = ginit; + dsave[(2 - (1)) + _dsave_offset] = gtest; + dsave[(3 - (1)) + _dsave_offset] = gx; + dsave[(4 - (1)) + _dsave_offset] = gy; + dsave[(5 - (1)) + _dsave_offset] = finit; + dsave[(6 - (1)) + _dsave_offset] = fx; + dsave[(7 - (1)) + _dsave_offset] = fy; + dsave[(8 - (1)) + _dsave_offset] = stx; + dsave[(9 - (1)) + _dsave_offset] = sty; + dsave[(10 - (1)) + _dsave_offset] = stmin; + dsave[(11 - (1)) + _dsave_offset] = stmax; + dsave[(12 - (1)) + _dsave_offset] = width; + dsave[(13 - (1)) + _dsave_offset] = width1; + } + + // + // c====================== The end of dcsrch ============================== + // + // c ********** + // c + // c Subroutine dcstep + // c + // c This subroutine computes a safeguarded step for a search + // c procedure and updates an interval that contains a step that + // c satisfies a sufficient decrease and a curvature condition. + // c + // c The parameter stx contains the step with the least function + // c value. If brackt is set to .true. then a minimizer has + // c been bracketed in an interval with endpoints stx and sty. + // c The parameter stp contains the current step. + // c The subroutine assumes that if brackt is set to .true. then + // c + // c min(stx,sty) < stp < max(stx,sty), + // c + // c and that the derivative at stx is negative in the direction + // c of the step. + // c + // c The subroutine statement is + // c + // c subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt, + // c stpmin,stpmax) + // c + // c where + // c + // c stx is a double precision variable. + // c On entry stx is the best step obtained so far and is an + // c endpoint of the interval that contains the minimizer. + // c On exit stx is the updated best step. + // c + // c fx is a double precision variable. + // c On entry fx is the function at stx. + // c On exit fx is the function at stx. + // c + // c dx is a double precision variable. + // c On entry dx is the derivative of the function at + // c stx. The derivative must be negative in the direction of + // c the step, that is, dx and stp - stx must have opposite + // c signs. + // c On exit dx is the derivative of the function at stx. + // c + // c sty is a double precision variable. + // c On entry sty is the second endpoint of the interval that + // c contains the minimizer. + // c On exit sty is the updated endpoint of the interval that + // c contains the minimizer. + // c + // c fy is a double precision variable. + // c On entry fy is the function at sty. + // c On exit fy is the function at sty. + // c + // c dy is a double precision variable. + // c On entry dy is the derivative of the function at sty. + // c On exit dy is the derivative of the function at the exit sty. + + // c + // c stp is a double precision variable. + // c On entry stp is the current step. If brackt is set to .true. + // c then on input stp must be between stx and sty. + // c On exit stp is a new trial step. + // c + // c fp is a double precision variable. + // c On entry fp is the function at stp + // c On exit fp is unchanged. + // c + // c dp is a double precision variable. + // c On entry dp is the the derivative of the function at stp. + // c On exit dp is unchanged. + // c + // c brackt is an logical variable. + // c On entry brackt specifies if a minimizer has been bracketed. + // c Initially brackt must be set to .false. + // c On exit brackt specifies if a minimizer has been bracketed. + // c When a minimizer is bracketed brackt is set to .true. + // c + // c stpmin is a double precision variable. + // c On entry stpmin is a lower bound for the step. + // c On exit stpmin is unchanged. + // c + // c stpmax is a double precision variable. + // c On entry stpmax is an upper bound for the step. + // c On exit stpmax is unchanged. + // c + // c MINPACK-1 Project. June 1983 + // c Argonne National Laboratory. + // c Jorge J. More' and David J. Thuente. + // c + // c MINPACK-2 Project. October 1993. + // c Argonne National Laboratory and University of Minnesota. + // c Brett M. Averick and Jorge J. More'. + // c + // c ********** + // + // + internal static void dcstep(ref double stx, ref double fx, ref double dx, + ref double sty, ref double fy, ref double dy, ref double stp, + double fp, double dp, ref bool brackt, double stpmin, double stpmax) + { + double gamma = 0.0d; + double p = 0.0d; + double q = 0.0d; + double r = 0.0d; + double s = 0.0d; + double sgnd = 0.0d; + double stpc = 0.0d; + double stpf = 0.0d; + double stpq = 0.0d; + double theta = 0.0d; + sgnd = (dp * ((dx / System.Math.Abs(dx)))); + + // c First case: A higher function value. The minimum is bracketed. + // c If the cubic step is closer to stx than the quadratic step, the + // c cubic step is taken, otherwise the average of the cubic and + // c quadratic steps is taken. + + if ((fp > fx)) + { + theta = ((((3.0 * ((fx - fp))) / ((stp - stx))) + dx) + dp); + + s = System.Math.Max(System.Math.Abs(theta), System.Math.Max(System.Math.Abs(dx), System.Math.Abs(dp))); + + gamma = (s * System.Math.Sqrt(((System.Math.Pow(((theta / s)), 2)) - (((dx / s)) * ((dp / s)))))); + if ((stp < stx)) + { + gamma = (-(gamma)); + } + p = (((gamma - dx)) + theta); + q = (((((gamma - dx)) + gamma)) + dp); + r = (p / q); + stpc = (stx + (r * ((stp - stx)))); + stpq = (stx + (((((dx / (((((fx - fp)) / ((stp - stx))) + dx)))) / 2.0)) * ((stp - stx)))); + + if ((System.Math.Abs((stpc - stx)) < System.Math.Abs((stpq - stx)))) + { + stpf = stpc; + } + else + { + stpf = (stpc + (((stpq - stpc)) / 2.0)); + } + + brackt = true; + + // + // c Second case: A lower function value and derivatives of opposite + // c sign. The minimum is bracketed. If the cubic step is farther from + // c stp than the secant step, the cubic step is taken, otherwise the + + // c secant step is taken. + // + } + else if ((sgnd < 0.0)) + { + theta = ((((3.0 * ((fx - fp))) / ((stp - stx))) + dx) + dp); + s = System.Math.Max(System.Math.Abs(theta), System.Math.Max(System.Math.Abs(dx), System.Math.Abs(dp))); + gamma = (s * System.Math.Sqrt(((System.Math.Pow(((theta / s)), 2)) - (((dx / s)) * ((dp / s)))))); + + if ((stp > stx)) + { + gamma = (-(gamma)); + } + + p = (((gamma - dp)) + theta); + q = (((((gamma - dp)) + gamma)) + dx); + r = (p / q); + stpc = (stp + (r * ((stx - stp)))); + stpq = (stp + (((dp / ((dp - dx)))) * ((stx - stp)))); + + if ((System.Math.Abs((stpc - stp)) > System.Math.Abs((stpq - stp)))) + { + stpf = stpc; + } + else + { + stpf = stpq; + } + + brackt = true; + // + // c Third case: A lower function value, derivatives of the same sign, + + // c and the magnitude of the derivative decreases. + // + } + else if ((System.Math.Abs(dp) < System.Math.Abs(dx))) + { + // + // c The cubic step is computed only if the cubic tends to infinity + // c in the direction of the step or if the minimum of the cubic + // c is beyond stp. Otherwise the cubic step is defined to be the + // c secant step. + // + theta = ((((3.0 * ((fx - fp))) / ((stp - stx))) + dx) + dp); + s = System.Math.Max(System.Math.Abs(theta), System.Math.Max(System.Math.Abs(dx), System.Math.Abs(dp))); + + // + // c The case gamma = 0 only arises if the cubic does not tend + // c to infinity in the direction of the step. + // + gamma = (s * System.Math.Sqrt(System.Math.Max(0.0, + ((System.Math.Pow(((theta / s)), 2)) - (((dx / s)) * ((dp / s))))))); + + if ((stp > stx)) + { + gamma = (-(gamma)); + } + + p = (((gamma - dp)) + theta); + q = (((gamma + ((dx - dp)))) + gamma); + r = (p / q); + + if (((r < 0.0) && (gamma != 0.0))) + { + stpc = (stp + (r * ((stx - stp)))); + } + else if ((stp > stx)) + { + stpc = stpmax; + } + else + { + stpc = stpmin; + } + + stpq = (stp + (((dp / ((dp - dx)))) * ((stx - stp)))); + + if (brackt) + { + // + // c A minimizer has been bracketed. If the cubic step is + // c closer to stp than the secant step, the cubic step is + // c taken, otherwise the secant step is taken. + // + if ((System.Math.Abs((stpc - stp)) < System.Math.Abs((stpq - stp)))) + { + stpf = stpc; + } + else + { + stpf = stpq; + } + if ((stp > stx)) + { + stpf = System.Math.Min((stp + + (0.6600000000000000310862446895043831318617 * ((sty - stp)))), stpf); + } + else + { + stpf = System.Math.Max((stp + + (0.6600000000000000310862446895043831318617 * ((sty - stp)))), stpf); + } + } + else + { + // + // c A minimizer has not been bracketed. If the cubic step is + // c farther from stp than the secant step, the cubic step is + // c taken, otherwise the secant step is taken. + // + if ((System.Math.Abs((stpc - stp)) > System.Math.Abs((stpq - stp)))) + { + stpf = stpc; + } + else + { + stpf = stpq; + } + stpf = System.Math.Min(stpmax, stpf); + stpf = System.Math.Max(stpmin, stpf); + } + + // + // c Fourth case: A lower function value, derivatives of the same sign, + // c and the magnitude of the derivative does not decrease. If the + // c minimum is not bracketed, the step is either stpmin or stpmax, + // c otherwise the cubic step is taken. + // + } + else + { + if (brackt) + { + theta = ((((3.0 * ((fp - fy))) / ((sty - stp))) + dy) + dp); + s = System.Math.Max(System.Math.Abs(theta), System.Math.Max(System.Math.Abs(dy), System.Math.Abs(dp))); + gamma = (s * System.Math.Sqrt(((System.Math.Pow(((theta / s)), 2)) - (((dy / s)) * ((dp / s)))))); + + if ((stp > sty)) + { + gamma = (-(gamma)); + } + + p = (((gamma - dp)) + theta); + q = (((((gamma - dp)) + gamma)) + dy); + r = (p / q); + stpc = (stp + (r * ((sty - stp)))); + stpf = stpc; + } + else if ((stp > stx)) + { + stpf = stpmax; + } + else + { + stpf = stpmin; + } + } + + // + // c Update the interval which contains a minimizer. + // + if ((fp > fx)) + { + sty = stp; + fy = fp; + dy = dp; + } + else + { + if ((sgnd < 0.0)) + { + sty = stx; + fy = fx; + dy = dx; + } + stx = stp; + fx = fp; + dx = dp; + } + + // + // c Compute the new step. + // + stp = stpf; + } + + // + // c======================= The end of cmprlb ============================= + // + // + // + // c ************ + // c + // c Subroutine errclb + // c + // c This subroutine checks the validity of the input data. + // c + // c + // c * * * + // c + // c NEOS, November 1994. (Latest revision June 1996.) + // c Optimization Technology Center. + // c Argonne National Laboratory and Northwestern University. + // c Written by + // c Ciyou Zhu + // c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. + // c + // c + // c ************ + // + // + // c Check the input arguments for errors. + // + private static void errclb(int n, int m, double factr, + double[] l, int _l_offset, double[] u, int _u_offset, + int[] nbd, int _nbd_offset, ref Task task, ref int info, ref int k) + { + + int i = 0; + if ((n <= 0)) + { + task = Task.Error; + } + if ((m <= 0)) + { + task = Task.Error; + } + if ((factr < 0.0)) + { + task = Task.Error; + } + + // + // c Check the validity of the arrays nbd(i), u(i), and l(i). + // + { + for (i = 1; i <= n; i++) + { + if (((nbd[(i - (1)) + _nbd_offset] < 0) || (nbd[(i - (1)) + _nbd_offset] > 3))) + { + // c return + task = Task.Error; + info = -6; + k = i; + } + if ((nbd[(i - (1)) + _nbd_offset] == 2)) + { + if ((l[(i - (1)) + _l_offset] > u[(i - (1)) + _u_offset])) + { + // c return + task = Task.Error; + info = -7; + k = i; + } + } + } + } + + } + + // + // c======================= The end of errclb ============================= + // + // + // + // c ************ + // c + // c Subroutine formk + // c + // c This subroutine forms the LEL^T factorization of the indefinite + // c + // c matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] + // c [L_a -R_z theta*S'AA'S ] + // c where E = [-I 0] + // c [ 0 I] + // c The matrix K can be shown to be equal to the matrix M^[-1]N + // c occurring in section 5.1 of [1], as well as to the matrix + // c Mbar^[-1] Nbar in section 5.3. + // c + // c n is an integer variable. + // c On entry n is the dimension of the problem. + // c On exit n is unchanged. + // c + // c nsub is an integer variable + // c On entry nsub is the number of subspace variables in free set. + // c On exit nsub is not changed. + // c + // c ind is an integer array of dimension nsub. + // c On entry ind specifies the indices of subspace variables. + // c On exit ind is unchanged. + // c + // c nenter is an integer variable. + // c On entry nenter is the number of variables entering the + // c free set. + // c On exit nenter is unchanged. + // c + // c ileave is an integer variable. + // c On entry indx2(ileave),...,indx2(n) are the variables leaving + // c the free set. + // c On exit ileave is unchanged. + // c + // c indx2 is an integer array of dimension n. + // c On entry indx2(1),...,indx2(nenter) are the variables entering + // c the free set, while indx2(ileave),...,indx2(n) are the + // c variables leaving the free set. + // c On exit indx2 is unchanged. + // c + // c iupdat is an integer variable. + // c On entry iupdat is the total number of BFGS updates made so far. + // c On exit iupdat is unchanged. + // c + // c updatd is a logical variable. + // c On entry 'updatd' is true if the L-BFGS matrix is updatd. + // c On exit 'updatd' is unchanged. + // c + // c wn is a double precision array of dimension 2m x 2m. + // c On entry wn is unspecified. + // c On exit the upper triangle of wn stores the LEL^T factorization + + // c of the 2*col x 2*col indefinite matrix + // c [-D -Y'ZZ'Y/theta L_a'-R_z' ] + // c [L_a -R_z theta*S'AA'S ] + // c + // c wn1 is a double precision array of dimension 2m x 2m. + // c On entry wn1 stores the lower triangular part of + // c [Y' ZZ'Y L_a'+R_z'] + // c [L_a+R_z S'AA'S ] + // c in the previous iteration. + // c On exit wn1 stores the corresponding updated matrices. + // c The purpose of wn1 is just to store these inner products + // c so they can be easily updated and inserted into wn. + // c + // c m is an integer variable. + // c On entry m is the maximum number of variable metric corrections + + // c used to define the limited memory matrix. + // c On exit m is unchanged. + // c + // c ws, wy, sy, and wtyy are double precision arrays; + // c theta is a double precision variable; + // c col is an integer variable; + // c head is an integer variable. + // c On entry they store the information defining the + // c limited memory BFGS matrix: + // c ws(n,m) stores S, a set of s-vectors; + // c wy(n,m) stores Y, a set of y-vectors; + // c sy(m,m) stores S'Y; + // c wtyy(m,m) stores the Cholesky factorization + // c of (theta*S'S+LD^(-1)L') + // c theta is the scaling factor specifying B_0 = theta I; + // c col is the number of variable metric corrections stored; + // c head is the location of the 1st s- (or y-) vector in S (or Y). + // c On exit they are unchanged. + // c + // c info is an integer variable. + // c On entry info is unspecified. + // c On exit info = 0 for normal return; + // c = -1 when the 1st Cholesky factorization failed; + // c = -2 when the 2st Cholesky factorization failed. + // c + // c Subprograms called: + // c + // c Linpack ... dcopy, dpofa, + // c + // c + // c References: + // c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited + // c memory algorithm for bound constrained optimization'', + // c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. + // c + // c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a + // c limited memory FORTRAN code for solving bound constrained + // c optimization problems'', Tech. Report, NAM-11, EECS Department, + + // c Northwestern University, 1994. + // c + // c (Postscript files of these papers are available via anonymous + // c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) + // c + // c * * * + // c + // c NEOS, November 1994. (Latest revision June 1996.) + // c Optimization Technology Center. + // c Argonne National Laboratory and Northwestern University. + // c Written by + // c Ciyou Zhu + // c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. + // c + // c + // c ************ + // + // + // c Form the lower triangular part of + // c WN1 = [Y' ZZ'Y L_a'+R_z'] + // c [L_a+R_z S'AA'S ] + // c where L_a is the strictly lower triangular part of S'AA'Y + // c R_z is the upper triangular part of S'ZZ'Y. + // + private static void formk(int n, int nsub, int[] ind, int _ind_offset, + int nenter, int ileave, int[] indx2, int _indx2_offset, + int iupdat, bool updatd, double[] wn, int _wn_offset, double[] wn1, int _wn1_offset, + int m, double[] ws, int _ws_offset, double[] wy, int _wy_offset, + double[] sy, int _sy_offset, double theta, int col, int head, + ref int info) + { + + int m2 = 0; + int ipntr = 0; + int jpntr = 0; + int iy = 0; + int is2 = 0; + int jy = 0; + int js = 0; + int is1 = 0; + int js1 = 0; + int k1 = 0; + int i = 0; + int k = 0; + int col2 = 0; + int pbegin = 0; + int pend = 0; + int dbegin = 0; + int dend = 0; + int upcl = 0; + double temp1 = 0.0d; + double temp2 = 0.0d; + double temp3 = 0.0d; + double temp4 = 0.0d; + + if (updatd) + { + if ((iupdat > m)) + { + // c shift old part of WN1. + { + for (jy = 1; jy <= (m - 1); jy++) + { + js = (m + jy); + dcopy((m - jy), wn1, ((jy + 1) - (1)) + ((jy + 1) - (1)) * ((2 * m)) + + _wn1_offset, 1, wn1, (jy - (1)) + (jy - (1)) * ((2 * m)) + _wn1_offset, 1); + dcopy((m - jy), wn1, ((js + 1) - (1)) + ((js + 1) - (1)) * ((2 * m)) + + _wn1_offset, 1, wn1, (js - (1)) + (js - (1)) * ((2 * m)) + _wn1_offset, 1); + dcopy((m - 1), wn1, ((m + 2) - (1)) + ((jy + 1) - (1)) * ((2 * m)) + + _wn1_offset, 1, wn1, ((m + 1) - (1)) + (jy - (1)) * ((2 * m)) + _wn1_offset, 1); + } + } + } + + // + // c put new rows in blocks (1,1), (2,1) and (2,2). + pbegin = 1; + pend = nsub; + dbegin = (nsub + 1); + dend = n; + iy = col; + is2 = (m + col); + ipntr = ((head + col) - 1); + if ((ipntr > m)) + { + ipntr = (ipntr - m); + } + jpntr = head; + { + for (jy = 1; jy <= col; jy++) + { + js = (m + jy); + temp1 = 0.0; + temp2 = 0.0; + temp3 = 0.0; + + // c compute element jy of row 'col' of Y'ZZ'Y + { + for (k = pbegin; k <= pend; k++) + { + k1 = ind[(k - (1)) + _ind_offset]; + temp1 = (temp1 + (wy[(k1 - (1)) + (ipntr - (1)) + * (n) + _wy_offset] * wy[(k1 - (1)) + (jpntr - (1)) * (n) + _wy_offset])); + } + } + // c compute elements jy of row 'col' of L_a and S'AA'S + { + for (k = dbegin; k <= dend; k++) + { + k1 = ind[(k - (1)) + _ind_offset]; + temp2 = (temp2 + (ws[(k1 - (1)) + (ipntr - (1)) + * (n) + _ws_offset] * ws[(k1 - (1)) + (jpntr - (1)) * (n) + _ws_offset])); + temp3 = (temp3 + (ws[(k1 - (1)) + (ipntr - (1)) + * (n) + _ws_offset] * wy[(k1 - (1)) + (jpntr - (1)) * (n) + _wy_offset])); + } + } + + wn1[(iy - (1)) + (jy - (1)) * ((2 * m)) + _wn1_offset] = temp1; + wn1[(is2 - (1)) + (js - (1)) * ((2 * m)) + _wn1_offset] = temp2; + wn1[(is2 - (1)) + (jy - (1)) * ((2 * m)) + _wn1_offset] = temp3; + jpntr = ((jpntr) % (m) + 1); + } + } + + // + // c put new column in block (2,1). + jy = col; + jpntr = ((head + col) - 1); + if ((jpntr > m)) + { + jpntr = (jpntr - m); + } + ipntr = head; + { + for (i = 1; i <= col; i++) + { + is2 = (m + i); + temp3 = 0.0; + // c compute element i of column 'col' of R_z + { + for (k = pbegin; k <= pend; k++) + { + k1 = ind[(k - (1)) + _ind_offset]; + temp3 = (temp3 + (ws[(k1 - (1)) + (ipntr - (1)) + * (n) + _ws_offset] * wy[(k1 - (1)) + (jpntr - (1)) * (n) + _wy_offset])); + } + } + ipntr = ((ipntr) % (m) + 1); + wn1[(is2 - (1)) + (jy - (1)) * ((2 * m)) + _wn1_offset] = temp3; + } + } + upcl = (col - 1); + } + else + { + upcl = col; + } + + // + // c modify the old parts in blocks (1,1) and (2,2) due to changes + // c in the set of free variables. + + ipntr = head; + { + for (iy = 1; iy <= upcl; iy++) + { + is2 = (m + iy); + jpntr = head; + { + for (jy = 1; jy <= iy; jy++) + { + js = (m + jy); + temp1 = 0.0; + temp2 = 0.0; + temp3 = 0.0; + temp4 = 0.0; + { + for (k = 1; k <= nenter; k++) + { + k1 = indx2[(k - (1)) + _indx2_offset]; + temp1 = (temp1 + (wy[(k1 - (1)) + (ipntr - (1)) + * (n) + _wy_offset] * wy[(k1 - (1)) + (jpntr - (1)) * (n) + _wy_offset])); + temp2 = (temp2 + (ws[(k1 - (1)) + (ipntr - (1)) + * (n) + _ws_offset] * ws[(k1 - (1)) + (jpntr - (1)) * (n) + _ws_offset])); + } + } + { + for (k = ileave; k <= n; k++) + { + k1 = indx2[(k - (1)) + _indx2_offset]; + temp3 = (temp3 + (wy[(k1 - (1)) + (ipntr + - (1)) * (n) + _wy_offset] * wy[(k1 - (1)) + (jpntr - (1)) * (n) + _wy_offset])); + temp4 = (temp4 + (ws[(k1 - (1)) + (ipntr + - (1)) * (n) + _ws_offset] * ws[(k1 - (1)) + (jpntr - (1)) * (n) + _ws_offset])); + } + } + + wn1[(iy - (1)) + (jy - (1)) * ((2 * m)) + _wn1_offset] = + ((wn1[(iy - (1)) + (jy - (1)) * ((2 * m)) + _wn1_offset] + temp1) - temp3); + wn1[(is2 - (1)) + (js - (1)) * ((2 * m)) + _wn1_offset] = + ((wn1[(is2 - (1)) + (js - (1)) * ((2 * m)) + _wn1_offset] - temp2) + temp4); + jpntr = ((jpntr) % (m) + 1); + } + } + + ipntr = ((ipntr) % (m) + 1); + } + } + + // + // c modify the old parts in block (2,1). + ipntr = head; + { + for (is2 = (m + 1); is2 <= (m + upcl); is2++) + { + jpntr = head; + { + for (jy = 1; jy <= upcl; jy++) + { + temp1 = 0.0; + temp3 = 0.0; + { + for (k = 1; k <= nenter; k++) + { + k1 = indx2[(k - (1)) + _indx2_offset]; + temp1 = (temp1 + (ws[(k1 - (1)) + (ipntr + - (1)) * (n) + _ws_offset] * wy[(k1 - (1)) + (jpntr - (1)) * (n) + _wy_offset])); + } + } + { + for (k = ileave; k <= n; k++) + { + k1 = indx2[(k - (1)) + _indx2_offset]; + temp3 = (temp3 + (ws[(k1 - (1)) + (ipntr - (1)) + * (n) + _ws_offset] * wy[(k1 - (1)) + (jpntr - (1)) * (n) + _wy_offset])); + } + } + + if ((is2 <= (jy + m))) + { + wn1[(is2 - (1)) + (jy - (1)) * ((2 * m)) + _wn1_offset] = + ((wn1[(is2 - (1)) + (jy - (1)) * ((2 * m)) + _wn1_offset] + temp1) - temp3); + } + else + { + wn1[(is2 - (1)) + (jy - (1)) * ((2 * m)) + _wn1_offset] = + ((wn1[(is2 - (1)) + (jy - (1)) * ((2 * m)) + _wn1_offset] - temp1) + temp3); + } + jpntr = ((jpntr) % (m) + 1); + } + } + ipntr = ((ipntr) % (m) + 1); + } + } + + // + // c Form the upper triangle of WN = [D+Y' ZZ'Y/theta -L_a'+R_z' ] + // c [-L_a +R_z S'AA'S*theta] + // + m2 = (2 * m); + { + for (iy = 1; iy <= col; iy++) + { + is2 = (col + iy); + is1 = (m + iy); + { + for (jy = 1; jy <= iy; jy++) + { + js = (col + jy); + js1 = (m + jy); + wn[(jy - (1)) + (iy - (1)) * ((2 * m)) + _wn_offset] = + (wn1[(iy - (1)) + (jy - (1)) * ((2 * m)) + _wn1_offset] / theta); + wn[(js - (1)) + (is2 - (1)) * ((2 * m)) + _wn_offset] = + (wn1[(is1 - (1)) + (js1 - (1)) * ((2 * m)) + _wn1_offset] * theta); + } + } + { + for (jy = 1; jy <= (iy - 1); jy++) + { + wn[(jy - (1)) + (is2 - (1)) * ((2 * m)) + _wn_offset] = + (-(wn1[(is1 - (1)) + (jy - (1)) * ((2 * m)) + _wn1_offset])); + } + } + { + for (jy = iy; jy <= col; jy++) + { + wn[(jy - (1)) + (is2 - (1)) * ((2 * m)) + _wn_offset] = + wn1[(is1 - (1)) + (jy - (1)) * ((2 * m)) + _wn1_offset]; + } + } + + wn[(iy - (1)) + (iy - (1)) * ((2 * m)) + _wn_offset] = (wn[(iy - (1)) + (iy - (1)) * ((2 * m)) + + _wn_offset] + sy[(iy - (1)) + (iy - (1)) * (m) + _sy_offset]); + } + } + + // + // c Form the upper triangle of WN= [ LL' L^-1(-L_a'+R_z')] + // c [(-L_a +R_z)L'^-1 S'AA'S*theta ] + // + // c first Cholesky factor (1,1) block of wn to get LL' + // c with L' stored in the upper triangle of wn. + dpofa(wn, _wn_offset, m2, col, ref info); + + // c then form L^-1(-L_a'+R_z') in the (1,2) block. + col2 = (2 * col); + { + for (js = (col + 1); js <= col2; js++) + { + dtrsl(wn, _wn_offset, m2, col, wn, (1 - (1)) + + (js - (1)) * ((2 * m)) + _wn_offset, 11, ref info); + } + } + + // + // c Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the + // c upper triangle of (2,2) block of wn. + // + // + { + for (is2 = (col + 1); is2 <= col2; is2++) + { + { + for (js = is2; js <= col2; js++) + { + wn[(is2 - (1)) + (js - (1)) * ((2 * m)) + _wn_offset] = + (wn[(is2 - (1)) + (js - (1)) * ((2 * m)) + _wn_offset] + + ddot(col, wn, (1 - (1)) + (is2 - (1)) * ((2 * m)) + + _wn_offset, 1, wn, (1 - (1)) + (js - (1)) * ((2 * m)) + _wn_offset, 1)); + } + } + } + } + + // + // c Cholesky factorization of (2,2) block of wn. + // + dpofa(wn, ((col + 1) - (1)) + ((col + 1) - (1)) + * ((2 * m)) + _wn_offset, m2, col, ref info); + } + + // + // c======================= The end of formk ============================== + // + // + // + // c ************ + // c + // c Subroutine formt + // c + // c This subroutine forms the upper half of the pos. def. and symm. + + // c T = theta*SS + L*D^(-1)*L', stores T in the upper triangle + // c of the array wt, and performs the Cholesky factorization of T + + // c to produce J*J', with J' stored in the upper triangle of wt. + // c + // c Subprograms called: + // c + // c Linpack ... dpofa. + // c + // c + // c * * * + // c + // c NEOS, November 1994. (Latest revision June 1996.) + // c Optimization Technology Center. + // c Argonne National Laboratory and Northwestern University. + // c Written by + // c Ciyou Zhu + // c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. + // c + // c + // c ************ + // + // + // + // c Form the upper half of T = theta*SS + L*D^(-1)*L', + // c store T in the upper triangle of the array wt. + // + private static void formt(int m, + double[] wt, int _wt_offset, + double[] sy, int _sy_offset, + double[] ss, int _ss_offset, + int col, double theta, ref int info) + { + + int i = 0; + int j = 0; + int k = 0; + int k1 = 0; + + double ddum = 0.0d; + { + for (j = 1; j <= col; j++) + { + wt[(1 - (1)) + (j - (1)) * (m) + _wt_offset] = (theta * ss[(1 - (1)) + + (j - (1)) * (m) + _ss_offset]); + } + } + + { + for (i = 2; i <= col; i++) + { + { + for (j = i; j <= col; j++) + { + k1 = (System.Math.Min(i, j) - 1); + + ddum = 0.0; + { + for (k = 1; k <= k1; k++) + { + ddum = (ddum + ((sy[(i - (1)) + (k - (1)) * (m) + + _sy_offset] * sy[(j - (1)) + (k - (1)) * (m) + + _sy_offset]) / sy[(k - (1)) + (k - (1)) * (m) + _sy_offset])); + } + } + + wt[(i - (1)) + (j - (1)) * (m) + _wt_offset] = (ddum + + (theta * ss[(i - (1)) + (j - (1)) * (m) + _ss_offset])); + } + } + } + } + + // + // c Cholesky factorize T to J*J' with + // c J' stored in the upper triangle of wt. + // + dpofa(wt, _wt_offset, m, col, ref info); + } + + // + // c======================= The end of formt ============================== + // + // + // + // c ************ + // c + // c Subroutine freev + // c + // c This subroutine counts the entering and leaving variables when + // c iter > 0, and finds the index set of free and active variables + // c at the GCP. + // c + // c cnstnd is a logical variable indicating whether bounds are present + // c + // c index is an integer array of dimension n + // c for i=1,...,nfree, index(i) are the indices of free variables + // c for i=nfree+1,...,n, index(i) are the indices of bound variables + // c On entry after the first iteration, index gives + // c the free variables at the previous iteration. + // c On exit it gives the free variables based on the determination + // c in cauchy using the array iwhere. + // c + // c indx2 is an integer array of dimension n + // c On entry indx2 is unspecified. + // c On exit with iter>0, indx2 indicates which variables + // c have changed status since the previous iteration. + // c For i= 1,...,nenter, indx2(i) have changed from bound to free. + // c For i= ileave+1,...,n, indx2(i) have changed from free to bound. + // c + // c + // c * * * + // c + // c NEOS, November 1994. (Latest revision June 1996.) + // c Optimization Technology Center. + // c Argonne National Laboratory and Northwestern University. + // c Written by + // c Ciyou Zhu + // c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. + // c + // c + // c ************ + // + // + private static void freev(int n, ref int nfree, + int[] index, int _index_offset, ref int nenter, ref int ileave, + int[] indx2, int _indx2_offset, int[] iwhere, int _iwhere_offset, + ref bool wrk, bool updatd, bool cnstnd, int iprint, int iter) + { + + int iact = 0; + int i = 0; + int k = 0; + + nenter = 0; + ileave = (n + 1); + + if (((iter > 0) && cnstnd)) + { + // c count the entering and leaving variables. + { + for (i = 1; i <= nfree; i++) + { + k = index[(i - (1)) + _index_offset]; + + // + // c write(6,*) ' k = index(i) ', k + // c write(6,*) ' index = ', i + // + if ((iwhere[(k - (1)) + _iwhere_offset] > 0)) + { + ileave = (ileave - 1); + indx2[(ileave - (1)) + _indx2_offset] = k; + + if ((iprint >= 100)) + { + // DISPLAY: "Variable " + k + " leaves the set of free variables" + } + } + } + } + { + for (i = (1 + nfree); i <= n; i++) + { + k = index[(i - (1)) + _index_offset]; + if ((iwhere[(k - (1)) + _iwhere_offset] <= 0)) + { + nenter = (nenter + 1); + indx2[(nenter - (1)) + _indx2_offset] = k; + + + if ((iprint >= 100)) + { + // DISPLAY: "Variable " + k + " enters the set of free variables" + } + } + } + } + + if ((iprint >= 99)) + { + // DISPLAY: ((n + 1) - ileave)) + " variables leave; " + // nenter + " variables enter" + } + } + + wrk = ((((ileave < (n + 1))) || ((nenter > 0))) || updatd); + + // + // c Find the index set of free and active variables at the GCP. + // + nfree = 0; + iact = (n + 1); + { + for (i = 1; i <= n; i++) + { + if ((iwhere[(i - (1)) + _iwhere_offset] <= 0)) + { + nfree = (nfree + 1); + index[(nfree - (1)) + _index_offset] = i; + } + else + { + iact = (iact - 1); + index[(iact - (1)) + _index_offset] = i; + } + } + } + + if ((iprint >= 99)) + { + // DISPLAY: nfree + " variables are free at GCP " + (iter + 1)) + } + } + + // + // c======================= The end of freev ============================== + // + // + // c ************ + // c + // c Subroutine hpsolb + // c + // c This subroutine sorts out the least element of t, and puts the + // c remaining elements of t in a heap. + // c + // c n is an integer variable. + // c On entry n is the dimension of the arrays t and iorder. + // c On exit n is unchanged. + // c + // c t is a double precision array of dimension n. + // c On entry t stores the elements to be sorted, + // c On exit t(n) stores the least elements of t, and t(1) to t(n-1) + + // c stores the remaining elements in the form of a heap. + // c + // c iorder is an integer array of dimension n. + // c On entry iorder(i) is the index of t(i). + // c On exit iorder(i) is still the index of t(i), but iorder may be + + // c permuted in accordance with t. + // c + // c iheap is an integer variable specifying the task. + // c On entry iheap should be set as follows: + // c iheap .eq. 0 if t(1) to t(n) is not in the form of a heap, + // c iheap .ne. 0 if otherwise. + // c On exit iheap is unchanged. + // c + // c + // c References: + // c Algorithm 232 of CACM (J. W. J. Williams): HEAPSORT. + // c + // c * * * + // c + // c NEOS, November 1994. (Latest revision June 1996.) + // c Optimization Technology Center. + // c Argonne National Laboratory and Northwestern University. + // c Written by + // c Ciyou Zhu + // c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. + // c + // c ************ + // + // + private static void hpsolb(int n, double[] t, int _t_offset, + int[] iorder, int _iorder_offset, int iheap) + { + + int i = 0; + int j = 0; + int k = 0; + int indxin = 0; + int indxou = 0; + double ddum = 0.0d; + double out2 = 0.0d; + + if ((iheap == 0)) + { + // + // c Rearrange the elements t(1) to t(n) to form a heap. + // + { + for (k = 2; k <= n; k++) + { + ddum = t[(k - (1)) + _t_offset]; + indxin = iorder[(k - (1)) + _iorder_offset]; + // + // c Add ddum to the heap. + i = k; + + L10: + if ((i > 1)) + { + j = (i / 2); + if ((ddum < t[(j - (1)) + _t_offset])) + { + t[(i - (1)) + _t_offset] = t[(j - (1)) + _t_offset]; + iorder[(i - (1)) + _iorder_offset] = iorder[(j - (1)) + _iorder_offset]; + i = j; + goto L10; + } + } + + t[(i - (1)) + _t_offset] = ddum; + iorder[(i - (1)) + _iorder_offset] = indxin; + } + } + } + // + // c Assign to 'out' the value of t(1), the least member of the heap, + // c and rearrange the remaining members to form a heap as + // c elements 1 to n-1 of t. + // + if ((n > 1)) + { + i = 1; + out2 = t[(1 - (1)) + _t_offset]; + indxou = iorder[(1 - (1)) + _iorder_offset]; + ddum = t[(n - (1)) + _t_offset]; + indxin = iorder[(n - (1)) + _iorder_offset]; + + // + // c Restore the heap + L30: + j = (i + i); + if ((j <= (n - 1))) + { + if ((t[((j + 1) - (1)) + _t_offset] < t[(j - (1)) + _t_offset])) + { + j = (j + 1); + } + if ((t[(j - (1)) + _t_offset] < ddum)) + { + t[(i - (1)) + _t_offset] = t[(j - (1)) + _t_offset]; + iorder[(i - (1)) + _iorder_offset] = iorder[(j - (1)) + _iorder_offset]; + i = j; + goto L30; + } + } + t[(i - (1)) + _t_offset] = ddum; + iorder[(i - (1)) + _iorder_offset] = indxin; + // + // c Put the least member in t(n). + // + t[(n - (1)) + _t_offset] = out2; + iorder[(n - (1)) + _iorder_offset] = indxou; + } + } + + // + // c====================== The end of hpsolb ============================== + // + // + // c ********** + // c + // c Subroutine lnsrlb + // c + // c This subroutine calls subroutine dcsrch from the Minpack2 library + + // c to perform the line search. Subroutine dscrch is safeguarded so + // c that all trial points lie within the feasible region. + // c + // c Subprograms called: + // c + // c Minpack2 Library ... dcsrch. + // c + // c Linpack ... dtrsl, + // c + // c + // c * * * + // c + // c NEOS, November 1994. (Latest revision June 1996.) + // c Optimization Technology Center. + // c Argonne National Laboratory and Northwestern University. + // c Written by + // c Ciyou Zhu + // c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. + // c + // c + // c ********** + // + // + private static void lnsrlb(int n, double[] l, int _l_offset, double[] u, int _u_offset, + int[] nbd, int _nbd_offset, double[] x, int _x_offset, double f, ref double fold, ref double gd, + ref double gdold, double[] g, int _g_offset, double[] d, int _d_offset, double[] r, int _r_offset, + double[] t, int _t_offset, double[] z, int _z_offset, ref double stp, + ref double dnorm, ref double dtd, ref double xstep, ref double stpmx, int iter, + ref int ifun, ref int iback, ref int nfgv, ref int info, ref Task task, + bool boxed, bool cnstnd, ref Task csave, int[] isave, int _isave_offset, + double[] dsave, int _dsave_offset) + { + + int i = 0; + double a1 = 0.0d; + double a2 = 0.0d; + + if (task == Task.FG_LN) + { + goto L556; + } + + // + dtd = ddot(n, d, _d_offset, 1, d, _d_offset, 1); + dnorm = System.Math.Sqrt(dtd); + // + // c Determine the maximum step length. + // + stpmx = 10000000000.0; + if (cnstnd) + { + if ((iter == 0)) + { + stpmx = 1.0; + } + else + { + { + for (i = 1; i <= n; i++) + { + a1 = d[(i - (1)) + _d_offset]; + if ((nbd[(i - (1)) + _nbd_offset] != 0)) + { + if (((a1 < 0.0) && (nbd[(i - (1)) + _nbd_offset] <= 2))) + { + a2 = (l[(i - (1)) + _l_offset] - x[(i - (1)) + _x_offset]); + if ((a2 >= 0.0)) + { + stpmx = 0.0; + } + else if (((a1 * stpmx) < a2)) + { + stpmx = (a2 / a1); + } + } + else if (((a1 > 0.0) && (nbd[(i - (1)) + _nbd_offset] >= 2))) + { + a2 = (u[(i - (1)) + _u_offset] - x[(i - (1)) + _x_offset]); + if ((a2 <= 0.0)) + { + stpmx = 0.0; + } + else if (((a1 * stpmx) > a2)) + { + stpmx = (a2 / a1); + } + } + } + } + } + } + } + + // + if (((iter == 0) && (!boxed))) + { + if (double.IsNaN(dnorm)) + stp = stpmx; + else + stp = System.Math.Min((1.0 / dnorm), stpmx); + } + else + { + stp = 1.0; + } + + // + dcopy(n, x, _x_offset, 1, t, _t_offset, 1); + dcopy(n, g, _g_offset, 1, r, _r_offset, 1); + + fold = f; + ifun = 0; + iback = 0; + csave = Task.Start; + + L556: + gd = BFGS.ddot(n, g, _g_offset, 1, d, _d_offset, 1); + + if ((ifun == 0)) + { + gdold = gd; + if ((gd >= 0.0)) + { + // the directional derivative >=0. + // Line search is impossible. + + // DISPLAY: " ascent direction in projection gd = " + gd + info = -4; + return; + } + } + + dcsrch(f, gd, ref stp, + 0.001000000000000000020816681711721685132943, + 0.9000000000000000222044604925031308084726, + 0.1000000000000000055511151231257827021182, 0.0, + stpmx, ref csave, isave, _isave_offset, dsave, _dsave_offset); + + // + xstep = (stp * dnorm); + if (csave != Task.Convergence && csave != Task.Warning) + { + task = Task.FG_LN; + ifun = (ifun + 1); + nfgv = (nfgv + 1); + iback = (ifun - 1); + if ((stp == 1.0)) + { + dcopy(n, z, _z_offset, 1, x, _x_offset, 1); + } + else + { + { + for (i = 1; i <= n; i++) + { + x[(i - (1)) + _x_offset] = ((stp * d[(i - (1)) + _d_offset]) + t[(i - (1)) + _t_offset]); + } + } + } + } + else + { + task = Task.New_X; + } + } + + // + // c======================= The end of setulb ============================= + // + // c-jlm-jn + // + // + // c ************ + // c + // c Subroutine mainlb + // c + // c This subroutine solves bound constrained optimization problems by + + // c using the compact formula of the limited memory BFGS updates. + // c + // c n is an integer variable. + // c On entry n is the number of variables. + // c On exit n is unchanged. + // c + // c m is an integer variable. + // c On entry m is the maximum number of variable metric + // c corrections allowed in the limited memory matrix. + // c On exit m is unchanged. + // c + // c x is a double precision array of dimension n. + // c On entry x is an approximation to the solution. + // c On exit x is the current approximation. + // c + // c l is a double precision array of dimension n. + // c On entry l is the lower bound of x. + // c On exit l is unchanged. + // c + // c u is a double precision array of dimension n. + // c On entry u is the upper bound of x. + // c On exit u is unchanged. + // c + // c nbd is an integer array of dimension n. + // c On entry nbd represents the type of bounds imposed on the + // c variables, and must be specified as follows: + // c nbd(i)=0 if x(i) is unbounded, + // c 1 if x(i) has only a lower bound, + // c 2 if x(i) has both lower and upper bounds, + // c 3 if x(i) has only an upper bound. + // c On exit nbd is unchanged. + // c + // c f is a double precision variable. + // c On first entry f is unspecified. + // c On final exit f is the value of the function at x. + // c + // c g is a double precision array of dimension n. + // c On first entry g is unspecified. + // c On final exit g is the value of the gradient at x. + // c + // c factr is a double precision variable. + // c On entry factr >= 0 is specified by the user. The iteration + // c will stop when + // c + // c (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch + // c + // c where epsmch is the machine precision, which is automatically + + // c generated by the code. + // c On exit factr is unchanged. + // c + // c pgtol is a double precision variable. + // c On entry pgtol >= 0 is specified by the user. The iteration + // c will stop when + // c + // c max{|proj g_i | i = 1, ..., n} <= pgtol + // c + // c where pg_i is the ith component of the projected gradient. + // c On exit pgtol is unchanged. + // c + // c ws, wy, sy, and wt are double precision working arrays used to + // c store the following information defining the limited memory + // c BFGS matrix: + // c ws, of dimension n x m, stores S, the matrix of s-vectors; + // c wy, of dimension n x m, stores Y, the matrix of y-vectors; + // c sy, of dimension m x m, stores S'Y; + // c ss, of dimension m x m, stores S'S; + // c yy, of dimension m x m, stores Y'Y; + // c wt, of dimension m x m, stores the Cholesky factorization + // c of (theta*S'S+LD^(-1)L'); see eq. + // c (2.26) in [3]. + // c + // c wn is a double precision working array of dimension 2m x 2m + // c used to store the LEL^T factorization of the indefinite matrix + // c K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] + // c [L_a -R_z theta*S'AA'S ] + // c + // c where E = [-I 0] + // c [ 0 I] + // c + // c snd is a double precision working array of dimension 2m x 2m + // c used to store the lower triangular part of + // c N = [Y' ZZ'Y L_a'+R_z'] + // c [L_a +R_z S'AA'S ] + // c + // c z(n),r(n),d(n),t(n), xp(n),wa(8*m) are double precision working ar + // c z is used at different times to store the Cauchy point and + // c the Newton point. + // c xp is used to safeguard the projected Newton direction + // c + // c sg(m),sgo(m),yg(m),ygo(m) are double precision working arrays. + // c + // c index is an integer working array of dimension n. + // c In subroutine freev, index is used to store the free and fixed + // c variables at the Generalized Cauchy Point (GCP). + // c + // c iwhere is an integer working array of dimension n used to record + // c the status of the vector x for GCP computation. + // c iwhere(i)=0 or -3 if x(i) is free and has bounds, + // c 1 if x(i) is fixed at l(i), and l(i) .ne. u(i) + // c 2 if x(i) is fixed at u(i), and u(i) .ne. l(i) + // c 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i) + // c -1 if x(i) is always free, i.e., no bounds on it. + // c + // c indx2 is an integer working array of dimension n. + // c Within subroutine cauchy, indx2 corresponds to the array iorder. + // c In subroutine freev, a list of variables entering and leaving + // c the free set is stored in indx2, and it is passed on to + // c subroutine formk with this information. + // c + // c task is a working string of characters of length 60 indicating + // c the current job when entering and leaving this subroutine. + // c + // c iprint is an INTEGER variable that must be set by the user. + // c It controls the frequency and type of output generated: + // c iprint<0 no output is generated; + // c iprint=0 print only one line at the last iteration; + // c 0100 print details of every iteration including x and g; + // c When iprint > 0, the file iterate.dat will be created to + // c summarize the iteration. + // c + // c csave is a working string of characters of length 60. + // c + // c lsave is a logical working array of dimension 4. + // c + // c isave is an integer working array of dimension 23. + // c + // c dsave is a double precision working array of dimension 29. + // c + // c + // c Subprograms called + // c + // c L-BFGS-B Library ... cauchy, subsm, lnsrlb, formk, + // c + // c errclb, prn1lb, prn2lb, prn3lb, active, projgr, + // c + // c freev, cmprlb, matupd, formt. + // c + // c Minpack2 Library ... timer + // c + // c Linpack Library ... dcopy, + // c + // c + // c References: + // c + // c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited + // c memory algorithm for bound constrained optimization'', + // c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. + // c + // c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN + // c Subroutines for Large Scale Bound Constrained Optimization'' + // c Tech. Report, NAM-11, EECS Department, Northwestern University, + + // c 1994. + // c + // c [3] R. Byrd, J. Nocedal and R. Schnabel "Representations of + // c Quasi-Newton Matrices and their use in Limited Memory Methods'', + // c Mathematical Programming 63 (1994), no. 4, pp. 129-156. + // c + // c (Postscript files of these papers are available via anonymous + // c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) + // c + // c * * * + // c + // c NEOS, November 1994. (Latest revision June 1996.) + // c Optimization Technology Center. + // c Argonne National Laboratory and Northwestern University. + // c Written by + // c Ciyou Zhu + // c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. + // c + // c + // c ************ + // + // + private static void mainlb(int n, + int m, + double[] x, int _x_offset, + double[] l, int _l_offset, + double[] u, int _u_offset, + int[] nbd, int _nbd_offset, + ref double f, + double[] g, int _g_offset, + double factr, + double pgtol, + double[] ws, int _ws_offset, + double[] wy, int _wy_offset, + double[] sy, int _sy_offset, + double[] ss, int _ss_offset, + double[] wt, int _wt_offset, + double[] wn, int _wn_offset, + double[] snd, int _snd_offset, + double[] z, int _z_offset, + double[] r, int _r_offset, + double[] d, int _d_offset, + double[] t, int _t_offset, + double[] xp, int _xp_offset, + double[] wa, int _wa_offset, + int[] index, int _index_offset, + int[] iwhere, int _iwhere_offset, + int[] indx2, int _indx2_offset, + ref Task task, + int iprint, + ref Task csave, + bool[] lsave, int _lsave_offset, + int[] isave, int _isave_offset, + double[] dsave, int _dsave_offset) + { + + bool prjctd = false; + bool cnstnd = false; + bool boxed = false; + bool updatd = false; + bool wrk = false; + int i = 0; + int k = 0; + int nintol = 0; + int itfile = 0; + int iback = 0; + int nskip = 0; + int head = 0; + int col = 0; + int iter = 0; + int itail = 0; + int iupdat = 0; + int nseg = 0; + int nfgv = 0; + int info = 0; + int ifun = 0; + int iword = 0; + int nfree = 0; + int nact = 0; + int ileave = 0; + int nenter = 0; + double theta = 0.0d; + double fold = 0.0d; + double dr = 0.0d; + double rr = 0.0d; + double tol = 0.0d; + double xstep = 0.0d; + double sbgnrm = 0.0d; + double ddum = 0.0d; + double dnorm = 0.0d; + double dtd = 0.0d; + double epsmch = 0.0d; + double cpu1 = 0.0d; + double cpu2 = 0.0d; + double cachyt = 0.0d; + double sbtime = 0.0d; + double lnscht = 0.0d; + double time1 = 0.0d; + // double time2 = 0.0d; + double gd = 0.0d; + double gdold = 0.0d; + double stp = 0.0d; + double stpmx = 0.0d; + // double time = 0.0d; + // float epsilon = 0.0f; + + if (task == Task.Start) + { + // + // epsmch = (double)(Epsilon.epsilon(1.0)); + epsmch = 1.11022302462516E-16; + + // + //Timer.timer(time1); + // + // c Initialize counters and scalars when task='START'. + // + // c for the limited memory BFGS matrices: + col = 0; + head = 1; + theta = 1.0; + iupdat = 0; + updatd = false; + iback = 0; + itail = 0; + iword = 0; + nact = 0; + ileave = 0; + nenter = 0; + fold = 0.0; + dnorm = 0.0; + cpu1 = 0.0; + gd = 0.0; + stpmx = 0.0; + sbgnrm = 0.0; + stp = 0.0; + gdold = 0.0; + dtd = 0.0; + // + // c for operation counts: + iter = 0; + nfgv = 0; + nseg = 0; + nintol = 0; + nskip = 0; + nfree = n; + ifun = 0; + // c for stopping tolerance: + tol = (factr * epsmch); + // + // c for measuring running time: + cachyt = (double)(0); + sbtime = (double)(0); + lnscht = (double)(0); + + // + // c 'info' records the termination information. + info = 0; + // + itfile = 8; + if ((iprint >= 1)) + { + // c open a summary file 'iterate.dat' + ; // WARNING: Unimplemented statement in Fortran source. + } + // + // c Check the input arguments for errors. + // + errclb(n, m, factr, l, _l_offset, u, _u_offset, + nbd, _nbd_offset, ref task, ref info, ref k); + + if (task == Task.Error) + { + /* + Prn3lb.prn3lb(n, x, _x_offset, f, task, iprint, + info, itfile, iter, nfgv, nintol, nskip, nact, sbgnrm, 0.0, nseg, word, + iback, stp, xstep, k, cachyt, sbtime, lnscht); + */ + return; + } + + // + // Prn1lb.prn1lb(n, m, l, _l_offset, u, _u_offset, x, _x_offset, iprint, itfile, epsmch); + + // + // c Initialize iwhere & project x onto the feasible set. + // + active(n, l, _l_offset, u, _u_offset, nbd, _nbd_offset, x, + _x_offset, iwhere, _iwhere_offset, iprint, ref prjctd, ref cnstnd, ref boxed); + // + // c The end of the initialization. + // + } + else + { + // c restore local variables. + // + prjctd = lsave[(1 - (1)) + _lsave_offset]; + cnstnd = lsave[(2 - (1)) + _lsave_offset]; + boxed = lsave[(3 - (1)) + _lsave_offset]; + updatd = lsave[(4 - (1)) + _lsave_offset]; + // + nintol = isave[(1 - (1)) + _isave_offset]; + itfile = isave[(3 - (1)) + _isave_offset]; + iback = isave[(4 - (1)) + _isave_offset]; + nskip = isave[(5 - (1)) + _isave_offset]; + head = isave[(6 - (1)) + _isave_offset]; + col = isave[(7 - (1)) + _isave_offset]; + itail = isave[(8 - (1)) + _isave_offset]; + iter = isave[(9 - (1)) + _isave_offset]; + iupdat = isave[(10 - (1)) + _isave_offset]; + nseg = isave[(12 - (1)) + _isave_offset]; + nfgv = isave[(13 - (1)) + _isave_offset]; + info = isave[(14 - (1)) + _isave_offset]; + ifun = isave[(15 - (1)) + _isave_offset]; + iword = isave[(16 - (1)) + _isave_offset]; + nfree = isave[(17 - (1)) + _isave_offset]; + nact = isave[(18 - (1)) + _isave_offset]; + ileave = isave[(19 - (1)) + _isave_offset]; + nenter = isave[(20 - (1)) + _isave_offset]; + // + theta = dsave[(1 - (1)) + _dsave_offset]; + fold = dsave[(2 - (1)) + _dsave_offset]; + tol = dsave[(3 - (1)) + _dsave_offset]; + dnorm = dsave[(4 - (1)) + _dsave_offset]; + epsmch = dsave[(5 - (1)) + _dsave_offset]; + cpu1 = dsave[(6 - (1)) + _dsave_offset]; + cachyt = dsave[(7 - (1)) + _dsave_offset]; + sbtime = dsave[(8 - (1)) + _dsave_offset]; + lnscht = dsave[(9 - (1)) + _dsave_offset]; + time1 = dsave[(10 - (1)) + _dsave_offset]; + gd = dsave[(11 - (1)) + _dsave_offset]; + stpmx = dsave[(12 - (1)) + _dsave_offset]; + sbgnrm = dsave[(13 - (1)) + _dsave_offset]; + stp = dsave[(14 - (1)) + _dsave_offset]; + gdold = dsave[(15 - (1)) + _dsave_offset]; + dtd = dsave[(16 - (1)) + _dsave_offset]; + // + // c After returning from the driver go to the point where execution + // c is to resume. + // + if (task == Task.FG_LN) + { + goto L666; + } + + if (task == Task.New_X) + { + goto L777; + } + + if (task == Task.FG_ST) + { + goto L111; + } + } + + // + // c Compute f0 and g0. + // + task = Task.FG_ST; + + // c return to the driver to calculate f and g; reenter at 111. + goto L1000; + + L111: + + nfgv = 1; + + // + // c Compute the infinity norm of the (-) projected gradient. + // + projgr(n, l, _l_offset, u, _u_offset, nbd, _nbd_offset, + x, _x_offset, g, _g_offset, ref sbgnrm); + + if ((sbgnrm <= pgtol)) + { + // terminate the algorithm. + task = Task.Convergence; + goto L999; + } + // + // c ----------------- the beginning of the loop -------------------------- + // + L222: + iword = -1; + //Compute the Generalized Cauchy Point (GCP). + cauchy(n, x, _x_offset, l, _l_offset, u, _u_offset, nbd, _nbd_offset, + g, _g_offset, indx2, _indx2_offset, iwhere, _iwhere_offset, t, _t_offset, + d, _d_offset, z, _z_offset, m, wy, _wy_offset, ws, _ws_offset, sy, _sy_offset, + wt, _wt_offset, theta, col, head, wa, (1 - (1)) + _wa_offset, wa, + (((2 * m) + 1) - (1)) + _wa_offset, wa, (((4 * m) + 1) - (1)) + _wa_offset, wa, + (((6 * m) + 1) - (1)) + _wa_offset, ref nseg, iprint, sbgnrm, ref info, epsmch); + + cachyt = cachyt + cpu2 - cpu1; + nintol = nintol + nseg; + // Count the entering and leaving variables for iter > 0; + // find the index set of free and active variables at the GCP. + freev(n, ref nfree, index, _index_offset, ref nenter, ref ileave, indx2, _indx2_offset, + iwhere, _iwhere_offset, ref wrk, updatd, cnstnd, iprint, iter); + nact = (n - nfree); + // If there are no free variables or B=theta*I, + // then skip the subspace minimization. + // + if (nfree == 0 || col == 0) { + goto L555; + } + // c Subspace minimization. + // c Form the LEL^T factorization of the indefinite + // c matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] + // c [L_a -R_z theta*S'AA'S ] + // c where E = [-I 0] + // c [ 0 I] + if (wrk) { + formk(n, nfree, index, _index_offset, nenter, + ileave, indx2, _indx2_offset, iupdat, updatd, wn, _wn_offset, + snd, _snd_offset, m, ws, _ws_offset, wy, _wy_offset, sy, _sy_offset, + theta, col, head, ref info); + } + // compute r=-Z'B(xcp-xk)-Z'g (using wa(2m+1)=W'(xcp-x) from 'cauchy') + cmprlb(n, m, x, _x_offset, g, _g_offset, ws, _ws_offset, wy, _wy_offset, + sy, _sy_offset, wt, _wt_offset, z, _z_offset, r, _r_offset, wa, _wa_offset, + index, _index_offset, theta, col, head, nfree, cnstnd, ref info); + // c-jlm-jn call the direct method. + subsm(n, m, nfree, index, _index_offset, l, _l_offset, u, _u_offset, + nbd, _nbd_offset, z, _z_offset, r, _r_offset, xp, _xp_offset, ws, _ws_offset, + wy, _wy_offset, theta, x, _x_offset, g, _g_offset, col, head, + ref iword, wa, _wa_offset, wn, _wn_offset, iprint, ref info); + + L555: + // c Line search and optimality tests. + // c Generate the search direction d:=z-x. + for (i = 1; i <= n; i++) { + d[(i - (1)) + _d_offset] = (z[(i - (1)) + _z_offset] - x[(i - (1)) + _x_offset]); + } + + L666: + lnsrlb(n, l, _l_offset, u, _u_offset, nbd, _nbd_offset, x, _x_offset, + f, ref fold, ref gd, ref gdold, g, _g_offset, d, _d_offset, r, _r_offset, t, _t_offset, + z, _z_offset, ref stp, ref dnorm, ref dtd, ref xstep, ref stpmx, iter, ref ifun, + ref iback, ref nfgv, ref info, ref task, boxed, cnstnd, ref csave, isave, + (22 - (1)) + _isave_offset, dsave, (17 - (1)) + _dsave_offset); + + if (iback >= 20) + { + // restore the previous iterate. + dcopy(n, t, _t_offset, 1, x, _x_offset, 1); + dcopy(n, r, _r_offset, 1, g, _g_offset, 1); + f = fold; + if ((col == 0)) + { + info = -9; + // restore the actual number of f and g evaluations etc. + nfgv = (nfgv - 1); + ifun = (ifun - 1); + iback = (iback - 1); + task = Task.Abnormal; + iter = (iter + 1); + goto L999; + } + else + { + nfgv = (nfgv - 1); + info = 0; + col = 0; + head = 1; + theta = 1.0; + iupdat = 0; + updatd = false; + task = Task.Restart_LN; + goto L222; + } + } + else if (task == Task.FG_LN) + { + // return to the driver for calculating f and g; reenter at 666. + goto L1000; + } + else + { + // calculate and print out the quantities related to the new X. + iter = (iter + 1); + // Compute the infinity norm of the projected (-)gradient. + projgr(n, l, _l_offset, u, _u_offset, nbd, _nbd_offset, x, _x_offset, g, _g_offset, ref sbgnrm); + goto L1000; + } + + L777: + // c Test for termination. + if ((sbgnrm <= pgtol)) + { + // terminate the algorithm. + task = Task.Convergence; + goto L999; + } + + ddum = System.Math.Max(System.Math.Abs(fold), System.Math.Max(System.Math.Abs(f), 1.0)); + + if ((((fold - f)) <= (tol * ddum))) + { + // terminate the algorithm. + task = Task.Convergence; + + if ((iback >= 10)) + { + info = -5; + } + + // i.e., to issue a warning if iback>10 in the line search. + goto L999; + } + + // + // Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's. + // + for (i = 1; i <= n; i++) + { + r[(i - (1)) + _r_offset] = (g[(i - (1)) + _g_offset] - r[(i - (1)) + _r_offset]); + } + + rr = ddot(n, r, _r_offset, 1, r, _r_offset, 1); + + if ((stp == 1.0)) + { + dr = (gd - gdold); + ddum = (-(gdold)); + } + else + { + dr = (((gd - gdold)) * stp); + dscal(n, stp, d, _d_offset, 1); + ddum = (-((gdold * stp))); + } + + if ((dr <= (epsmch * ddum))) + { + // skip the L-BFGS update. + nskip = (nskip + 1); + updatd = false; + + + if ((iprint >= 1)) + { + // DISPLAY: dr, ddum + // DISPLAY: ' ys=',1p,e10.3,' -gs=',1P,E10.3,' BFGS update SKIPPED'" + } + + goto L888; + } + + // + // cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + // c + // c Update the L-BFGS matrix. + // c + // cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + // + updatd = true; + iupdat = (iupdat + 1); + + // + // c Update matrices WS and WY and form the middle matrix in B. + // + matupd(n, m, ws, _ws_offset, wy, _wy_offset, sy, _sy_offset, + ss, _ss_offset, d, _d_offset, r, _r_offset, ref itail, iupdat, ref col, + ref head, ref theta, rr, dr, stp, dtd); + + // + // c Form the upper half of the pds T = theta*SS + L*D^(-1)*L'; + // c Store T in the upper triangular of the array wt; + // c Cholesky factorize T to J*J' with + // c J' stored in the upper triangular of wt. + // + formt(m, wt, _wt_offset, sy, _sy_offset, ss, + _ss_offset, col, theta, ref info); + + // + // Now the inverse of the middle matrix in B is + // + // [ D^(1/2) O ] [ -D^(1/2) D^(-1/2)*L' ] + // [ -L*D^(-1/2) J ] [ 0 J' ] + // + L888: + // + // c -------------------- the end of the loop ----------------------------- + // + goto L222; + + L999: + + + L1000: + // + // Save local variables. + // + lsave[(1 - (1)) + _lsave_offset] = prjctd; + lsave[(2 - (1)) + _lsave_offset] = cnstnd; + lsave[(3 - (1)) + _lsave_offset] = boxed; + lsave[(4 - (1)) + _lsave_offset] = updatd; + // + isave[(1 - (1)) + _isave_offset] = nintol; + isave[(3 - (1)) + _isave_offset] = itfile; + isave[(4 - (1)) + _isave_offset] = iback; + isave[(5 - (1)) + _isave_offset] = nskip; + isave[(6 - (1)) + _isave_offset] = head; + isave[(7 - (1)) + _isave_offset] = col; + isave[(8 - (1)) + _isave_offset] = itail; + isave[(9 - (1)) + _isave_offset] = iter; + isave[(10 - (1)) + _isave_offset] = iupdat; + isave[(12 - (1)) + _isave_offset] = nseg; + isave[(13 - (1)) + _isave_offset] = nfgv; + isave[(14 - (1)) + _isave_offset] = info; + isave[(15 - (1)) + _isave_offset] = ifun; + isave[(16 - (1)) + _isave_offset] = iword; + isave[(17 - (1)) + _isave_offset] = nfree; + isave[(18 - (1)) + _isave_offset] = nact; + isave[(19 - (1)) + _isave_offset] = ileave; + isave[(20 - (1)) + _isave_offset] = nenter; + // + dsave[(1 - (1)) + _dsave_offset] = theta; + dsave[(2 - (1)) + _dsave_offset] = fold; + dsave[(3 - (1)) + _dsave_offset] = tol; + dsave[(4 - (1)) + _dsave_offset] = dnorm; + dsave[(5 - (1)) + _dsave_offset] = epsmch; + dsave[(6 - (1)) + _dsave_offset] = cpu1; + dsave[(7 - (1)) + _dsave_offset] = cachyt; + dsave[(8 - (1)) + _dsave_offset] = sbtime; + dsave[(9 - (1)) + _dsave_offset] = lnscht; + dsave[(10 - (1)) + _dsave_offset] = time1; + dsave[(11 - (1)) + _dsave_offset] = gd; + dsave[(12 - (1)) + _dsave_offset] = stpmx; + dsave[(13 - (1)) + _dsave_offset] = sbgnrm; + dsave[(14 - (1)) + _dsave_offset] = stp; + dsave[(15 - (1)) + _dsave_offset] = gdold; + dsave[(16 - (1)) + _dsave_offset] = dtd; + + return; + } + + // + // c======================= The end of lnsrlb ============================= + // + // + // + // c ************ + // c + // c Subroutine matupd + // c + // c This subroutine updates matrices WS and WY, and forms the + // c middle matrix in B. + // c + // c Subprograms called: + // c + // c Linpack ... dcopy, + // c + // c + // c * * * + // c + // c NEOS, November 1994. (Latest revision June 1996.) + // c Optimization Technology Center. + // c Argonne National Laboratory and Northwestern University. + // c Written by + // c Ciyou Zhu + // c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. + // c + // c + // c ************ + // + // + // c Set pointers for matrices WS and WY. + // + internal static void matupd(int n, + int m, double[] ws, int _ws_offset, double[] wy, int _wy_offset, + double[] sy, int _sy_offset, double[] ss, int _ss_offset, double[] d, int _d_offset, + double[] r, int _r_offset, ref int itail, int iupdat, ref int col, + ref int head, ref double theta, double rr, double dr, double stp, double dtd) + { + + int j = 0; + int pointr = 0; + + if ((iupdat <= m)) + { + col = iupdat; + itail = ((((head + iupdat) - 2)) % (m) + 1); + } + else + { + itail = ((itail) % (m) + 1); + head = ((head) % (m) + 1); + } + + // + // c Update matrices WS and WY. + // + dcopy(n, d, _d_offset, 1, ws, (1 - (1)) + (itail - (1)) * (n) + _ws_offset, 1); + dcopy(n, r, _r_offset, 1, wy, (1 - (1)) + (itail - (1)) * (n) + _wy_offset, 1); + + // + // c Set theta=yy/ys. + // + theta = (rr / dr); + // + // c Form the middle matrix in B. + // + // c update the upper triangle of SS, + // c and the lower triangle of SY: + + if ((iupdat > m)) + { + // c move old information + { + for (j = 1; j <= (col - 1); j++) + { + dcopy(j, ss, (2 - (1)) + ((j + 1) - (1)) * (m) + + _ss_offset, 1, ss, (1 - (1)) + (j - (1)) * (m) + _ss_offset, 1); + dcopy((col - j), sy, ((j + 1) - (1)) + ((j + 1) + - (1)) * (m) + _sy_offset, 1, sy, (j - (1)) + (j - (1)) * (m) + _sy_offset, 1); + } + } + } + // c add new information: the last row of SY + // c and the last column of SS: + pointr = head; + { + for (j = 1; j <= (col - 1); j++) + { + sy[(col - (1)) + (j - (1)) * (m) + _sy_offset] = + ddot(n, + d, _d_offset, 1, + wy, (1 - (1)) + (pointr - (1)) * (n) + _wy_offset, 1); + + ss[(j - (1)) + (col - (1)) * (m) + _ss_offset] = + ddot(n, + ws, (1 - (1)) + (pointr - (1)) * (n) + _ws_offset, 1, + d, _d_offset, 1); + + pointr = ((pointr) % (m) + 1); + } + } + if ((stp == 1.0)) + { + ss[(col - (1)) + (col - (1)) * (m) + _ss_offset] = dtd; + } + else + { + ss[(col - (1)) + (col - (1)) * (m) + _ss_offset] = ((stp * stp) * dtd); + } + + sy[(col - (1)) + (col - (1)) * (m) + _sy_offset] = dr; + + return; + } + + // + // c======================= The end of prn3lb ============================= + // + // + // + // c ************ + // c + // c Subroutine projgr + // c + // c This subroutine computes the infinity norm of the projected + // c gradient. + // c + // c + // c * * * + // c + // c NEOS, November 1994. (Latest revision June 1996.) + // c Optimization Technology Center. + // c Argonne National Laboratory and Northwestern University. + // c Written by + // c Ciyou Zhu + // c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. + // c + // c + // c ************ + // + // + internal static void projgr(int n, + double[] l, int _l_offset, + double[] u, int _u_offset, + int[] nbd, int _nbd_offset, + double[] x, int _x_offset, + double[] g, int _g_offset, + ref double sbgnrm) + { + + int i = 0; + double gi = 0.0d; + sbgnrm = 0.0; + { + for (i = 1; i <= n; i++) + { + gi = g[(i - (1)) + _g_offset]; + if ((nbd[(i - (1)) + _nbd_offset] != 0)) + { + if ((gi < 0.0)) + { + if ((nbd[(i - (1)) + _nbd_offset] >= 2)) + { + gi = System.Math.Max(((x[(i - (1)) + _x_offset] - u[(i - (1)) + _u_offset])), gi); + } + } + else + { + if ((nbd[(i - (1)) + _nbd_offset] <= 2)) + { + gi = System.Math.Min(((x[(i - (1)) + _x_offset] - l[(i - (1)) + _l_offset])), gi); + } + } + } + sbgnrm = System.Math.Max(sbgnrm, System.Math.Abs(gi)); + } + } + } + + // c + // c L-BFGS-B is released under the �New BSD License� (aka �Modified + // c or �3-clause license�) + // c Please read attached file License.txt + // c + // c=========== L-BFGS-B (version 3.0. April 25, 2011 ================= + // c + // c This is a modified version of L-BFGS-B. Minor changes in the updat + // c code appear preceded by a line comment as follows + // c + // c c-jlm-jn + // c + // c Major changes are described in the accompanying paper: + // c + // c Jorge Nocedal and Jose Luis Morales, Remark on "Algorithm 778: + // c L-BFGS-B: Fortran Subroutines for Large-Scale Bound Constraine + // c Optimization" (2011). To appear in ACM Transactions on + // c Mathematical Software, + // c + // c The paper describes an improvement and a correction to Algorithm 7 + // c It is shown that the performance of the algorithm can be improved + // c significantly by making a relatively simple modication to the subs + // c minimization phase. The correction concerns an error caused by the + // c of routine dpmeps to estimate machine precision. + // c + // c The total work space **wa** required by the new version is + // c + // c 2*m*n + 11m*m + 5*n + 8*m + // c + // c the old version required + // c + // c 2*m*n + 12m*m + 4*n + 12*m + // c + // c + // c J. Nocedal Department of Electrical Engineering and + // c Computer Science. + // c Northwestern University. Evanston, IL. USA + // c + // c + // c J.L Morales Departamento de Matematicas, + // c Instituto Tecnologico Autonomo de Mexico + // c Mexico D.F. Mexico. + // c + // c March 2011 + // c + // c======================================================================= + // + // c + // c-jlm-jn + // + // + // c ************ + // c + // c Subroutine setulb + // c + // c This subroutine partitions the working arrays wa and iwa, and + // c then uses the limited memory BFGS method to solve the bound + // c constrained optimization problem by calling mainlb. + // c (The direct method will be used in the subspace minimization.) + // c + // c n is an integer variable. + // c On entry n is the dimension of the problem. + // c On exit n is unchanged. + // c + // c m is an integer variable. + // c On entry m is the maximum number of variable metric corrections + // c + // c used to define the limited memory matrix. + // c On exit m is unchanged. + // c + // c x is a double precision array of dimension n. + // c On entry x is an approximation to the solution. + // c On exit x is the current approximation. + // c + // c l is a double precision array of dimension n. + // c On entry l is the lower bound on x. + // c On exit l is unchanged. + // c + // c u is a double precision array of dimension n. + // c On entry u is the upper bound on x. + // c On exit u is unchanged. + // c + // c nbd is an integer array of dimension n. + // c On entry nbd represents the type of bounds imposed on the + // c variables, and must be specified as follows: + // c nbd(i)=0 if x(i) is unbounded, + // c 1 if x(i) has only a lower bound, + // c 2 if x(i) has both lower and upper bounds, and + // c 3 if x(i) has only an upper bound. + // c On exit nbd is unchanged. + // c + // c f is a double precision variable. + // c On first entry f is unspecified. + // c On final exit f is the value of the function at x. + // c + // c g is a double precision array of dimension n. + // c On first entry g is unspecified. + // c On final exit g is the value of the gradient at x. + // c + // c factr is a double precision variable. + // c On entry factr >= 0 is specified by the user. The iteration + // c will stop when + // c + // c (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch + // c + // c where epsmch is the machine precision, which is automatically + // c generated by the code. Typical values for factr: 1.d+12 for + // c low accuracy; 1.d+7 for moderate accuracy; 1.d+1 for extremely + // c high accuracy. + // c On exit factr is unchanged. + // c + // c pgtol is a double precision variable. + // c On entry pgtol >= 0 is specified by the user. The iteration + // c will stop when + // c + // c max{|proj g_i | i = 1, ..., n} <= pgtol + // c + // c where pg_i is the ith component of the projected gradient. + + // c On exit pgtol is unchanged. + // c + // c wa is a double precision working array of length + // c (2mmax + 5)nmax + 12mmax^2 + 12mmax. + // c + // c iwa is an integer working array of length 3nmax. + // c + // c task is a working string of characters of length 60 indicating + // c the current job when entering and quitting this subroutine. + // c + // c iprint is an integer variable that must be set by the user. + // c It controls the frequency and type of output generated: + // c iprint<0 no output is generated; + // c iprint=0 print only one line at the last iteration; + // c 0100 print details of every iteration including x and g; + // c When iprint > 0, the file iterate.dat will be created to + // c summarize the iteration. + // c + // c csave is a working string of characters of length 60. + // c + // c lsave is a logical working array of dimension 4. + // c On exit with 'task' = NEW_X, the following information is + // c available: + // c If lsave(1) = .true. then the initial X has been replaced by + // c its projection in the feasible set + // c If lsave(2) = .true. then the problem is constrained; + // c If lsave(3) = .true. then each variable has upper and lower + // c bounds; + // c + // c isave is an integer working array of dimension 44. + // c On exit with 'task' = NEW_X, the following information is + // c available: + // c isave(22) = the total number of intervals explored in the + // c search of Cauchy points; + // c isave(26) = the total number of skipped BFGS updates before + // c the current iteration; + // c isave(30) = the number of current iteration; + // c isave(31) = the total number of BFGS updates prior the current + // c iteration; + // c isave(33) = the number of intervals explored in the search of + // c Cauchy point in the current iteration; + // c isave(34) = the total number of function and gradient + // c evaluations; + // c isave(36) = the number of function value or gradient + // c evaluations in the current iteration; + // c if isave(37) = 0 then the subspace argmin is within the box; + // c if isave(37) = 1 then the subspace argmin is beyond the box; + // c isave(38) = the number of free variables in the current + // c iteration; + // c isave(39) = the number of active constraints in the current + // c iteration; + // c n + 1 - isave(40) = the number of variables leaving the set of + // c active constraints in the current iteration; + // c isave(41) = the number of variables entering the set of active + // c constraints in the current iteration. + // c + // c dsave is a double precision working array of dimension 29. + // c On exit with 'task' = NEW_X, the following information is + // c available: + // c dsave(1) = current 'theta' in the BFGS matrix; + // c dsave(2) = f(x) in the previous iteration; + // c dsave(3) = factr*epsmch; + // c dsave(4) = 2-norm of the line search direction vector; + // c dsave(5) = the machine precision epsmch generated by the code; + // c dsave(7) = the accumulated time spent on searching for + // c Cauchy points; + // c dsave(8) = the accumulated time spent on + // c subspace minimization; + // c dsave(9) = the accumulated time spent on line search; + // c dsave(11) = the slope of the line search function at + // c the current point of line search; + // c dsave(12) = the maximum relative step length imposed in + // c line search; + // c dsave(13) = the infinity norm of the projected gradient; + // c dsave(14) = the relative step length in the line search; + // c dsave(15) = the slope of the line search function at + // c the starting point of the line search; + // c dsave(16) = the square of the 2-norm of the line search + // c direction vector. + // c + // c Subprograms called: + // c + // c L-BFGS-B Library ... mainlb. + // c + // c + // c References: + // c + // c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited + // c memory algorithm for bound constrained optimization'', + // c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. + // c + // c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a + // c limited memory FORTRAN code for solving bound constrained + // c optimization problems'', Tech. Report, NAM-11, EECS Department, + // c Northwestern University, 1994. + // c + // c (Postscript files of these papers are available via anonymous + // c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) + // c + // c * * * + // c + // c NEOS, November 1994. (Latest revision June 1996.) + // c Optimization Technology Center. + // c Argonne National Laboratory and Northwestern University. + // c Written by + // c Ciyou Zhu + // c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. + // c + // c + // c ************ + // + internal static void setulb(int n, + int m, + double[] x, int _x_offset, + double[] l, int _l_offset, + double[] u, int _u_offset, + int[] nbd, int _nbd_offset, + ref double f, + double[] g, int _g_offset, + double factr, + double pgtol, + double[] wa, int _wa_offset, + int[] iwa, int _iwa_offset, + ref Task task, + int iprint, + ref Task csave, + bool[] lsave, int _lsave_offset, + int[] isave, int _isave_offset, + double[] dsave, int _dsave_offset) + { + + int lws = 0; + int lr = 0; + int lz = 0; + int lt = 0; + int ld = 0; + int lxp = 0; + int lwa = 0; + int lwy = 0; + int lsy = 0; + int lss = 0; + int lwt = 0; + int lwn = 0; + int lsnd = 0; + + if (task == Task.Start) + { + isave[(1 - (1)) + _isave_offset] = (m * n); + isave[(2 - (1)) + _isave_offset] = ((int)System.Math.Pow(m, 2)); + isave[(3 - (1)) + _isave_offset] = (4 * ((int)System.Math.Pow(m, 2))); + isave[(4 - (1)) + _isave_offset] = 1; + isave[(5 - (1)) + _isave_offset] = (isave[(4 - (1)) + _isave_offset] + isave[(1 - (1)) + _isave_offset]); + isave[(6 - (1)) + _isave_offset] = (isave[(5 - (1)) + _isave_offset] + isave[(1 - (1)) + _isave_offset]); + isave[(7 - (1)) + _isave_offset] = (isave[(6 - (1)) + _isave_offset] + isave[(2 - (1)) + _isave_offset]); + isave[(8 - (1)) + _isave_offset] = (isave[(7 - (1)) + _isave_offset] + isave[(2 - (1)) + _isave_offset]); + isave[(9 - (1)) + _isave_offset] = (isave[(8 - (1)) + _isave_offset] + isave[(2 - (1)) + _isave_offset]); + isave[(10 - (1)) + _isave_offset] = (isave[(9 - (1)) + _isave_offset] + isave[(3 - (1)) + _isave_offset]); + isave[(11 - (1)) + _isave_offset] = (isave[(10 - (1)) + _isave_offset] + isave[(3 - (1)) + _isave_offset]); + isave[(12 - (1)) + _isave_offset] = (isave[(11 - (1)) + _isave_offset] + n); + isave[(13 - (1)) + _isave_offset] = (isave[(12 - (1)) + _isave_offset] + n); + isave[(14 - (1)) + _isave_offset] = (isave[(13 - (1)) + _isave_offset] + n); + isave[(15 - (1)) + _isave_offset] = (isave[(14 - (1)) + _isave_offset] + n); + isave[(16 - (1)) + _isave_offset] = (isave[(15 - (1)) + _isave_offset] + n); + } + lws = isave[(4 - (1)) + _isave_offset]; + lwy = isave[(5 - (1)) + _isave_offset]; + lsy = isave[(6 - (1)) + _isave_offset]; + lss = isave[(7 - (1)) + _isave_offset]; + lwt = isave[(8 - (1)) + _isave_offset]; + lwn = isave[(9 - (1)) + _isave_offset]; + lsnd = isave[(10 - (1)) + _isave_offset]; + lz = isave[(11 - (1)) + _isave_offset]; + lr = isave[(12 - (1)) + _isave_offset]; + ld = isave[(13 - (1)) + _isave_offset]; + lt = isave[(14 - (1)) + _isave_offset]; + lxp = isave[(15 - (1)) + _isave_offset]; + lwa = isave[(16 - (1)) + _isave_offset]; + // + mainlb(n, m, x, _x_offset, l, _l_offset, u, _u_offset, nbd, + _nbd_offset, ref f, g, _g_offset, factr, pgtol, wa, (lws - (1)) + _wa_offset, + wa, (lwy - (1)) + _wa_offset, wa, (lsy - (1)) + _wa_offset, wa, + (lss - (1)) + _wa_offset, wa, (lwt - (1)) + _wa_offset, wa, + (lwn - (1)) + _wa_offset, wa, (lsnd - (1)) + _wa_offset, wa, + (lz - (1)) + _wa_offset, wa, (lr - (1)) + _wa_offset, wa, + (ld - (1)) + _wa_offset, wa, (lt - (1)) + _wa_offset, wa, + (lxp - (1)) + _wa_offset, wa, (lwa - (1)) + _wa_offset, iwa, + (1 - (1)) + _iwa_offset, iwa, ((n + 1) - (1)) + _iwa_offset, + iwa, (((2 * n) + 1) - (1)) + _iwa_offset, ref task, iprint, ref csave, + lsave, _lsave_offset, isave, (22 - (1)) + _isave_offset, dsave, + _dsave_offset); + } + + // + // c======================= The end of projgr ============================= + // + // + // c ****************************************************************** + // c + // c This routine contains the major changes in the updated version. + // c The changes are described in the accompanying paper + // c + // c Jose Luis Morales, Jorge Nocedal + // c "Remark On Algorithm 788: L-BFGS-B: Fortran Subroutines for Large + // c Bound Constrained Optimization". Decemmber 27, 2010. + // c + // c J.L. Morales Departamento de Matematicas, + // c Instituto Tecnologico Autonomo de Mexico + // c Mexico D.F. + // c + // c J, Nocedal Department of Electrical Engineering and + // c Computer Science. + // c Northwestern University. Evanston, IL. USA + // c + // c January 17, 2011 + // c + // c ***************************************************************** + // c + // c + // c Subroutine subsm + // c + // c Given xcp, l, u, r, an index set that specifies + // c the active set at xcp, and an l-BFGS matrix B + // c (in terms of WY, WS, SY, WT, head, col, and theta), + // c this subroutine computes an approximate solution + // c of the subspace problem + // c + // c (P) min Q(x) = r'(x-xcp) + 1/2 (x-xcp)' B (x-xcp) + // c + // c subject to l<=x<=u + // c x_i=xcp_i for all i in A(xcp) + // c + // c along the subspace unconstrained Newton direction + // c + // c d = -(Z'BZ)^(-1) r. + // c + // c The formula for the Newton direction, given the L-BFGS matrix + // c and the Sherman-Morrison formula, is + // c + // c d = (1/theta)r + (1/theta*2) Z'WK^(-1)W'Z r. + // c + // c where + // c K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] + // c [L_a -R_z theta*S'AA'S ] + // c + // c Note that this procedure for computing d differs + // c from that described in [1]. One can show that the matrix K is + // c equal to the matrix M^[-1]N in that paper. + // c + // c n is an integer variable. + // c On entry n is the dimension of the problem. + // c On exit n is unchanged. + // c + // c m is an integer variable. + // c On entry m is the maximum number of variable metric corrections + + // c used to define the limited memory matrix. + // c On exit m is unchanged. + // c + // c nsub is an integer variable. + // c On entry nsub is the number of free variables. + // c On exit nsub is unchanged. + // c + // c ind is an integer array of dimension nsub. + // c On entry ind specifies the coordinate indices of free variables. + // c On exit ind is unchanged. + // c + // c l is a double precision array of dimension n. + // c On entry l is the lower bound of x. + // c On exit l is unchanged. + // c + // c u is a double precision array of dimension n. + // c On entry u is the upper bound of x. + // c On exit u is unchanged. + // c + // c nbd is a integer array of dimension n. + // c On entry nbd represents the type of bounds imposed on the + // c variables, and must be specified as follows: + // c nbd(i)=0 if x(i) is unbounded, + // c 1 if x(i) has only a lower bound, + // c 2 if x(i) has both lower and upper bounds, and + // c 3 if x(i) has only an upper bound. + // c On exit nbd is unchanged. + // c + // c x is a double precision array of dimension n. + // c On entry x specifies the Cauchy point xcp. + // c On exit x(i) is the minimizer of Q over the subspace of + // c free variables. + // c + // c d is a double precision array of dimension n. + // c On entry d is the reduced gradient of Q at xcp. + // c On exit d is the Newton direction of Q. + // c + // c xp is a double precision array of dimension n. + // c used to safeguard the projected Newton direction + // c + // c xx is a double precision array of dimension n + // c On entry it holds the current iterate + // c On output it is unchanged + // + // c gg is a double precision array of dimension n + // c On entry it holds the gradient at the current iterate + // c On output it is unchanged + // c + // c ws and wy are double precision arrays; + // c theta is a double precision variable; + // c col is an integer variable; + // c head is an integer variable. + // c On entry they store the information defining the + // c limited memory BFGS matrix: + // c ws(n,m) stores S, a set of s-vectors; + // c wy(n,m) stores Y, a set of y-vectors; + // c theta is the scaling factor specifying B_0 = theta I; + // c col is the number of variable metric corrections stored; + // c head is the location of the 1st s- (or y-) vector in S (or Y). + // c On exit they are unchanged. + // c + // c iword is an integer variable. + // c On entry iword is unspecified. + // c On exit iword specifies the status of the subspace solution. + // c iword = 0 if the solution is in the box, + // c 1 if some bound is encountered. + // c + // c wv is a double precision working array of dimension 2m. + // c + // c wn is a double precision array of dimension 2m x 2m. + // c On entry the upper triangle of wn stores the LEL^T factorization + // c of the indefinite matrix + // c + // c K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] + // c [L_a -R_z theta*S'AA'S ] + // c where E = [-I 0] + // c [ 0 I] + // c On exit wn is unchanged. + // c + // c iprint is an INTEGER variable that must be set by the user. + // c It controls the frequency and type of output generated: + // c iprint<0 no output is generated; + // c iprint=0 print only one line at the last iteration; + // c 0100 print details of every iteration including x and g; + // c When iprint > 0, the file iterate.dat will be created to + // c summarize the iteration. + // c + // c info is an integer variable. + // c On entry info is unspecified. + // c On exit info = 0 for normal return, + // c = nonzero for abnormal return + // c when the matrix K is ill-conditioned. + // c + // c Subprograms called: + // c + // c Linpack + // c + // c + // c References: + // c + // c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited + // c memory algorithm for bound constrained optimization'', + // c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. + // c + // c + // c + // c * * * + // c + // c NEOS, November 1994. (Latest revision June 1996.) + // c Optimization Technology Center. + // c Argonne National Laboratory and Northwestern University. + // c Written by + // c Ciyou Zhu + // c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. + // c + // c + // c ************ + // + // c + // + internal static void subsm(int n, + int m, int nsub, int[] ind, int _ind_offset, + double[] l, int _l_offset, double[] u, int _u_offset, + int[] nbd, int _nbd_offset, double[] x, int _x_offset, + double[] d, int _d_offset, double[] xp, int _xp_offset, + double[] ws, int _ws_offset, double[] wy, int _wy_offset, + double theta, double[] xx, int _xx_offset, double[] gg, int _gg_offset, + int col, int head, ref int iword, + double[] wv, int _wv_offset, double[] wn, int _wn_offset, int iprint, ref int info) + { + int pointr = 0; + int m2 = 0; + int col2 = 0; + int ibd = 0; + int jy = 0; + int js = 0; + int i = 0; + int j = 0; + int k = 0; + double alpha = 0.0d; + double xk = 0.0d; + double dk = 0.0d; + double temp1 = 0.0d; + double temp2 = 0.0d; + double dd_p = 0.0d; + + if ((nsub <= 0)) + { + return; + } + + // + // Compute wv = W'Zd. + // + pointr = head; + { + for (i = 1; i <= col; i++) + { + temp1 = 0.0; + temp2 = 0.0; + { + for (j = 1; j <= nsub; j++) + { + k = ind[(j - (1)) + _ind_offset]; + temp1 = (temp1 + (wy[(k - (1)) + (pointr - (1)) + * (n) + _wy_offset] * d[(j - (1)) + _d_offset])); + temp2 = (temp2 + (ws[(k - (1)) + (pointr - (1)) + * (n) + _ws_offset] * d[(j - (1)) + _d_offset])); + } + } + + wv[(i - (1)) + _wv_offset] = temp1; + wv[((col + i) - (1)) + _wv_offset] = (theta * temp2); + pointr = ((pointr) % (m) + 1); + } + } + + // + // Compute wv:=K^(-1)wv. + // + m2 = (2 * m); + col2 = (2 * col); + + dtrsl(wn, _wn_offset, m2, col2, wv, _wv_offset, 11, ref info); + + { + for (i = 1; i <= col; i++) + { + wv[(i - (1)) + _wv_offset] = (-(wv[(i - (1)) + _wv_offset])); + } + } + + dtrsl(wn, _wn_offset, m2, col2, wv, _wv_offset, 1, ref info); + + // + // Compute d = (1/theta)d + (1/theta**2)Z'W wv. + // + pointr = head; + { + for (jy = 1; jy <= col; jy++) + { + js = (col + jy); + { + for (i = 1; i <= nsub; i++) + { + k = ind[(i - (1)) + _ind_offset]; + d[(i - (1)) + _d_offset] = ((d[(i - (1)) + _d_offset] + + ((wy[(k - (1)) + (pointr - (1)) * (n) + _wy_offset] * wv[(jy - (1)) + + _wv_offset]) / theta)) + (ws[(k - (1)) + + (pointr - (1)) * (n) + _ws_offset] * wv[(js - (1)) + _wv_offset])); + } + } + + pointr = ((pointr) % (m) + 1); + } + } + + + dscal(nsub, (1.0 / theta), d, _d_offset, 1); + + // + // ---------------------------------------------------- + // Let us try the projection, d is the Newton direction + // + iword = 0; + + dcopy(n, x, _x_offset, 1, xp, _xp_offset, 1); + + // c + { + for (i = 1; i <= nsub; i++) + { + k = ind[(i - (1)) + _ind_offset]; + dk = d[(i - (1)) + _d_offset]; + xk = x[(k - (1)) + _x_offset]; + if ((nbd[(k - (1)) + _nbd_offset] != 0)) + { + if ((nbd[(k - (1)) + _nbd_offset] == 1)) + { + x[(k - (1)) + _x_offset] = System.Math.Max(l[(k - (1)) + _l_offset], (xk + dk)); + if ((x[(k - (1)) + _x_offset] == l[(k - (1)) + _l_offset])) + { + iword = 1; + } + } + else + { + if ((nbd[(k - (1)) + _nbd_offset] == 2)) + { + xk = System.Math.Max(l[(k - (1)) + _l_offset], (xk + dk)); + x[(k - (1)) + _x_offset] = System.Math.Min(u[(k - (1)) + _u_offset], xk); + if (((x[(k - (1)) + _x_offset] == l[(k - (1)) + + _l_offset]) || (x[(k - (1)) + _x_offset] == u[(k - (1)) + _u_offset]))) + { + iword = 1; + } + } + else + { + if ((nbd[(k - (1)) + _nbd_offset] == 3)) + { + x[(k - (1)) + _x_offset] = System.Math.Min(u[(k - (1)) + _u_offset], (xk + dk)); + if ((x[(k - (1)) + _x_offset] == u[(k - (1)) + _u_offset])) + { + iword = 1; + } + } + } + } + } + else + { + x[(k - (1)) + _x_offset] = (xk + dk); + } + } + } + + + if ((iword == 0)) + { + goto L911; + } + + // + // check sign of the directional derivative + // + dd_p = 0.0; + { + for (i = 1; i <= n; i++) + { + dd_p = (dd_p + (((x[(i - (1)) + _x_offset] + - xx[(i - (1)) + _xx_offset])) * gg[(i - (1)) + _gg_offset])); + } + } + + if ((dd_p > 0.0)) + { + dcopy(n, xp, _xp_offset, 1, x, _x_offset, 1); + + // DISPLAY: " Positive dir derivative in projection " + // DISPLAY: " Using the backtracking step " + } + else + { + goto L911; + } + + // + // ----------------------------------------------------------------- + // + alpha = 1.0; + temp1 = alpha; + ibd = 0; + { + for (i = 1; i <= nsub; i++) + { + k = ind[(i - (1)) + _ind_offset]; + dk = d[(i - (1)) + _d_offset]; + if ((nbd[(k - (1)) + _nbd_offset] != 0)) + { + if (((dk < 0.0) && (nbd[(k - (1)) + _nbd_offset] <= 2))) + { + temp2 = (l[(k - (1)) + _l_offset] - x[(k - (1)) + _x_offset]); + if ((temp2 >= 0.0)) + { + temp1 = 0.0; + } + else if (((dk * alpha) < temp2)) + { + temp1 = (temp2 / dk); + } + } + else if (((dk > 0.0) && (nbd[(k - (1)) + _nbd_offset] >= 2))) + { + temp2 = (u[(k - (1)) + _u_offset] - x[(k - (1)) + _x_offset]); + if ((temp2 <= 0.0)) + { + temp1 = 0.0; + } + else if (((dk * alpha) > temp2)) + { + temp1 = (temp2 / dk); + } + } + + if ((temp1 < alpha)) + { + alpha = temp1; + ibd = i; + } + } + } + } + + if ((alpha < 1.0)) + { + dk = d[(ibd - (1)) + _d_offset]; + k = ind[(ibd - (1)) + _ind_offset]; + if ((dk > 0.0)) + { + x[(k - (1)) + _x_offset] = u[(k - (1)) + _u_offset]; + d[(ibd - (1)) + _d_offset] = 0.0; + } + else if ((dk < 0.0)) + { + x[(k - (1)) + _x_offset] = l[(k - (1)) + _l_offset]; + d[(ibd - (1)) + _d_offset] = 0.0; + } + } + { + for (i = 1; i <= nsub; i++) + { + k = ind[(i - (1)) + _ind_offset]; + x[(k - (1)) + _x_offset] = (x[(k - (1)) + _x_offset] + (alpha * d[(i - (1)) + _d_offset])); + } + } + + + L911: + + if ((iprint >= 99)) + { + // DISPLAY: "----------------exit SUBSM --------------------" + } + + return; + } + + // c + // c L-BFGS-B is released under the �New BSD License� (aka �Modified + // c or �3-clause license�) + // c Please read attached file License.txt + // c + // c + // c dpofa factors a double precision symmetric positive definite + // c matrix. + // c + // c dpofa is usually called by dpoco, but it can be called + // c directly with a saving in time if rcond is not needed. + // c (time for dpoco) = (1 + 18/n)*(time for dpofa) . + // c + // c on entry + // c + // c a double precision(lda, n) + // c the symmetric matrix to be factored. only the + // c diagonal and upper triangle are used. + // c + // c lda integer + // c the leading dimension of the array a . + // c + // c n integer + // c the order of the matrix a . + // c + // c on return + // c + // c a an upper triangular matrix r so that a = trans(r)*r + + // c where trans(r) is the transpose. + // c the strict lower triangle is unaltered. + // c if info .ne. 0 , the factorization is not complete. + // c + // c info integer + // c = 0 for normal return. + // c = k signals an error condition. the leading minor + // c of order k is not positive definite. + // c + // c linpack. this version dated 08/14/78 . + // c cleve moler, university of new mexico, argonne national lab. + // c + // c subroutines and functions + // c + // c blas ddot + // c fortran sqrt + // c + // c internal variables + // c + // c begin block with ...exits to 40 + // c + // c + private static void dpofa(double[] a, int _a_offset, int lda, int n, ref int info) + { + double t = 0.0d; + double s = 0.0d; + int j = 0; + int jm1 = 0; + int k = 0; + { + for (j = 1; j <= n; j++) + { + info = j; + s = 0.0e0; + jm1 = (j - 1); + if ((jm1 < 1)) + { + goto L20; + } + + { + for (k = 1; k <= jm1; k++) + { + t = (a[(k - (1)) + (j - (1)) * (lda) + _a_offset] + - ddot((k - 1), a, (1 - (1)) + (k - (1)) * (lda) + + _a_offset, 1, a, (1 - (1)) + (j - (1)) * (lda) + _a_offset, 1)); + + t = (t / a[(k - (1)) + (k - (1)) * (lda) + _a_offset]); + a[(k - (1)) + (j - (1)) * (lda) + _a_offset] = t; + s = (s + (t * t)); + } + } + + L20: + s = (a[(j - (1)) + (j - (1)) * (lda) + _a_offset] - s); + // c ......exit + if ((s <= 0.0e0)) + { + goto L40; + } + + a[(j - (1)) + (j - (1)) * (lda) + _a_offset] = System.Math.Sqrt(s); + } + } + + info = 0; + + L40: + return; + } + + // * .. Scalar Arguments .. + // * .. + // * .. Array Arguments .. + // * .. + // * + // * Purpose + // * ======= + // ** + // * scales a vector by a constant. + // * uses unrolled loops for increment equal to one. + // * jack dongarra, linpack, 3/11/78. + // * modified 3/93 to return if incx .le. 0. + // * modified 12/3/93, array(1) declarations changed to array(*) + // * + // * + // * .. Local Scalars .. + // * .. + // * .. Intrinsic Functions .. + // * .. + private static void dscal(int n, double da, double[] dx, int _dx_offset, int incx) + { + + int i = 0; + int m = 0; + int mp1 = 0; + int nincx = 0; + + if (((n <= 0) || (incx <= 0))) + { + return; + } + + if ((incx == 1)) + { + goto L20; + } + + // * + // * code for increment not equal to 1 + // * + nincx = (n * incx); + { + int _i_inc = incx; + for (i = 1; (_i_inc < 0) ? i >= nincx : i <= nincx; i += _i_inc) + { + dx[(i - (1)) + _dx_offset] = (da * dx[(i - (1)) + _dx_offset]); + } + } + return; + // * + // * code for increment equal to 1 + // * + // * + // * clean-up loop + // * + L20: + m = (n) % (5); + + if ((m == 0)) + { + goto L40; + } + { + for (i = 1; i <= m; i++) + { + dx[(i - (1)) + _dx_offset] = (da * dx[(i - (1)) + _dx_offset]); + } + } + + if ((n < 5)) + { + return; + } + + L40: + mp1 = (m + 1); + { + int _i_inc = 5; + for (i = mp1; i <= n; i += _i_inc) + { + dx[(i - (1)) + _dx_offset] = (da * dx[(i - (1)) + _dx_offset]); + dx[((i + 1) - (1)) + _dx_offset] = (da * dx[((i + 1) - (1)) + _dx_offset]); + dx[((i + 2) - (1)) + _dx_offset] = (da * dx[((i + 2) - (1)) + _dx_offset]); + dx[((i + 3) - (1)) + _dx_offset] = (da * dx[((i + 3) - (1)) + _dx_offset]); + dx[((i + 4) - (1)) + _dx_offset] = (da * dx[((i + 4) - (1)) + _dx_offset]); + } + } + + return; + } + + // * .. Scalar Arguments .. + // * .. + // * .. Array Arguments .. + // * .. + // * + // * Purpose + // * ======= + // * + // * forms the dot product of two vectors. + // * uses unrolled loops for increments equal to one. + // * jack dongarra, linpack, 3/11/78. + // * modified 12/3/93, array(1) declarations changed to array(*) + // * + // * + // * .. Local Scalars .. + // * .. + // * .. Intrinsic Functions .. + // * .. + private static double ddot(int n, double[] dx, int _dx_offset, + int incx, double[] dy, int _dy_offset, int incy) + { + double dtemp = 0.0d; + int i = 0; + int ix = 0; + int iy = 0; + int m = 0; + int mp1 = 0; + double ddot = 0.0d; + ddot = 0.0e0; + dtemp = 0.0e0; + if ((n <= 0)) + { + return ddot; + } + + if (((incx == 1) && (incy == 1))) + { + goto L20; + } + // * + // * code for unequal increments or equal increments + // * not equal to 1 + // * + ix = 1; + iy = 1; + if ((incx < 0)) + { + ix = (((((-(n)) + 1)) * incx) + 1); + } + if ((incy < 0)) + { + iy = (((((-(n)) + 1)) * incy) + 1); + } + + { + for (i = 1; i <= n; i++) + { + dtemp = (dtemp + (dx[(ix - (1)) + _dx_offset] * dy[(iy - (1)) + _dy_offset])); + ix = (ix + incx); + iy = (iy + incy); + } + } + ddot = dtemp; + return ddot; + + // * + // * code for both increments equal to 1 + // * + // * + // * clean-up loop + // * + L20: + m = (n) % (5); + + if ((m == 0)) + { + goto L40; + } + + { + for (i = 1; i <= m; i++) + { + dtemp = (dtemp + (dx[(i - (1)) + _dx_offset] * dy[(i - (1)) + _dy_offset])); + } + } + + if ((n < 5)) + { + goto L60; + } + + L40: + mp1 = (m + 1); + { + int _i_inc = 5; + for (i = mp1; i <= n; i += _i_inc) + { + dtemp = (((((dtemp + (dx[(i - (1)) + _dx_offset] * dy[(i - (1)) + + _dy_offset])) + (dx[((i + 1) - (1)) + _dx_offset] * dy[((i + 1) - (1)) + + _dy_offset])) + (dx[((i + 2) - (1)) + _dx_offset] * dy[((i + 2) - (1)) + + _dy_offset])) + (dx[((i + 3) - (1)) + _dx_offset] * dy[((i + 3) - (1)) + + _dy_offset])) + (dx[((i + 4) - (1)) + _dx_offset] * dy[((i + 4) - (1)) + _dy_offset])); + } + } + L60: + ddot = dtemp; + return ddot; + } + + private static void dcopy(int n, double[] dx, int _dx_offset, int incx, + double[] dy, int _dy_offset, int incy) + { + + int i = 0; + int ix = 0; + int iy = 0; + int m = 0; + int mp1 = 0; + if ((n <= 0)) + { + return; + } + + if (((incx == 1) && (incy == 1))) + { + goto L20; + } + // c + // c code for unequal increments or equal increments + // c not equal to 1 + // c + ix = 1; + iy = 1; + if ((incx < 0)) + { + ix = (((((-(n)) + 1)) * incx) + 1); + } + if ((incy < 0)) + { + iy = (((((-(n)) + 1)) * incy) + 1); + } + { + for (i = 1; i <= n; i++) + { + dy[(iy - (1)) + _dy_offset] = dx[(ix - (1)) + _dx_offset]; + ix = (ix + incx); + iy = (iy + incy); + } + } + return; + // c + // c code for both increments equal to 1 + // c + // c + // c clean-up loop + // c + L20: + m = (n) % (7); + if ((m == 0)) + { + goto L40; + } + { + for (i = 1; i <= m; i++) + { + dy[(i - (1)) + _dy_offset] = dx[(i - (1)) + _dx_offset]; + } + } + if ((n < 7)) + { + return; + } + + L40: + mp1 = (m + 1); + { + int _i_inc = 7; + for (i = mp1; i <= n; i += _i_inc) + { + dy[(i - (1)) + _dy_offset] = dx[(i - (1)) + _dx_offset]; + dy[((i + 1) - (1)) + _dy_offset] = dx[((i + 1) - (1)) + _dx_offset]; + dy[((i + 2) - (1)) + _dy_offset] = dx[((i + 2) - (1)) + _dx_offset]; + dy[((i + 3) - (1)) + _dy_offset] = dx[((i + 3) - (1)) + _dx_offset]; + dy[((i + 4) - (1)) + _dy_offset] = dx[((i + 4) - (1)) + _dx_offset]; + dy[((i + 5) - (1)) + _dy_offset] = dx[((i + 5) - (1)) + _dx_offset]; + dy[((i + 6) - (1)) + _dy_offset] = dx[((i + 6) - (1)) + _dx_offset]; + } + } + } + + // * .. Scalar Arguments .. + // * .. + // * .. Array Arguments .. + // * .. + // * + // * Purpose + // * ======= + // * + // * DAXPY constant times a vector plus a vector. + // * uses unrolled loops for increments equal to one. + // * + // * Further Details + // * =============== + // * + // * jack dongarra, linpack, 3/11/78. + // * modified 12/3/93, array(1) declarations changed to array(*) + // * + // * ===================================================================== + // * + // * .. Local Scalars .. + // * .. + // * .. Intrinsic Functions .. + // * .. + private static void daxpy(int n, double da, + double[] dx, int _dx_offset, int incx, double[] dy, int _dy_offset, int incy) + { + + int i = 0; + int ix = 0; + int iy = 0; + int m = 0; + int mp1 = 0; + if ((n <= 0)) + { + return; + } + + if ((da == 0.0e0)) + { + return; + } + + if (((incx == 1) && (incy == 1))) + { + // * + // * code for both increments equal to 1 + // * + // * + // * clean-up loop + // * + m = (n) % (4); + if ((m != 0)) + { + { + for (i = 1; i <= m; i++) + { + dy[(i - (1)) + _dy_offset] = (dy[(i - (1)) + _dy_offset] + + (da * dx[(i - (1)) + _dx_offset])); + } + } + } + + if ((n < 4)) + { + return; + } + + mp1 = (m + 1); + { + int _i_inc = 4; + for (i = mp1; i <= n; i += _i_inc) + { + dy[(i - (1)) + _dy_offset] = (dy[(i - (1)) + _dy_offset] + + (da * dx[(i - (1)) + _dx_offset])); + dy[((i + 1) - (1)) + _dy_offset] = (dy[((i + 1) - (1)) + _dy_offset] + + (da * dx[((i + 1) - (1)) + _dx_offset])); + dy[((i + 2) - (1)) + _dy_offset] = (dy[((i + 2) - (1)) + _dy_offset] + + (da * dx[((i + 2) - (1)) + _dx_offset])); + dy[((i + 3) - (1)) + _dy_offset] = (dy[((i + 3) - (1)) + _dy_offset] + + (da * dx[((i + 3) - (1)) + _dx_offset])); + } + } + } + else + { + // * + // * code for unequal increments or equal increments + // * not equal to 1 + // * + ix = 1; + iy = 1; + if ((incx < 0)) + { + ix = (((((-(n)) + 1)) * incx) + 1); + } + if ((incy < 0)) + { + iy = (((((-(n)) + 1)) * incy) + 1); + } + { + for (i = 1; i <= n; i++) + { + dy[(iy - (1)) + _dy_offset] = (dy[(iy - (1)) + _dy_offset] + + (da * dx[(ix - (1)) + _dx_offset])); + ix = (ix + incx); + iy = (iy + incy); + } + } + } + return; + } + } + +} \ No newline at end of file diff --git a/Assets/External/BioIK/Helpers/Evolution.cs.meta b/Assets/External/BioIK/Helpers/Evolution.cs.meta new file mode 100644 index 00000000..c6e8690c --- /dev/null +++ b/Assets/External/BioIK/Helpers/Evolution.cs.meta @@ -0,0 +1,19 @@ +fileFormatVersion: 2 +guid: 293464419e12c4a878d56542a32530e6 +timeCreated: 1502296338 +licenseType: Store +MonoImporter: + serializedVersion: 2 + defaultReferences: [] + executionOrder: 0 + icon: {instanceID: 0} + userData: + assetBundleName: + assetBundleVariant: +AssetOrigin: + serializedVersion: 1 + productId: 67819 + packageName: Bio IK + packageVersion: 2.0d + assetPath: Assets/BioIK/Helpers/Evolution.cs + uploadId: 199125 diff --git a/Assets/External/BioIK/Helpers/Model.cs b/Assets/External/BioIK/Helpers/Model.cs new file mode 100644 index 00000000..27ba148a --- /dev/null +++ b/Assets/External/BioIK/Helpers/Model.cs @@ -0,0 +1,577 @@ +using UnityEngine; +using System.Collections.Generic; + +namespace BioIK { + public class Model { + //Reference to the character + private BioIK Character; + + //Reference to root + private BioSegment Root; + + //Offset to world + private double OPX, OPY, OPZ; //Offset rosition to world frame + private double ORX, ORY, ORZ, ORW; //Offset rotation to world frame + private double OSX, OSY, OSZ; //Offset scale to world Frame + + //Linked list of nodes in the model + public Node[] Nodes = new Node[0]; + + //Global pointers to the IK setup + public MotionPtr[] MotionPtrs = new MotionPtr[0]; + public ObjectivePtr[] ObjectivePtrs = new ObjectivePtr[0]; + + //Assigned Configuraton + private double[] Configuration; + private double[] Gradient; + private double[] Losses; + + //Simulated Configuration + private double[] PX,PY,PZ,RX,RY,RZ,RW; + private double[] SimulatedLosses; + + //Degree of Freedom + private int DoF; + + public Model(BioIK character) { + Character = character; + + //Set Root + Root = Character.FindSegment(Character.transform); + + //Create Root + AddNode(Root); + + //Build Model + BioObjective[] objectives = CollectObjectives(Root, new List()); + for(int i=0; i chain = Character.GetChain(Root, objectives[i].Segment); + for(int j=1; j objectives) { + for(int i=0; i reverseChain = new List(); + reverseChain.Add(Transform); + Node p = parent; + while(p != null) { + reverseChain.Add(p.Transform); + p = p.Parent; + } + reverseChain.Reverse(); + Chain = reverseChain.ToArray(); + } + + //Adds a child to this node + public void AddChild(Node child) { + System.Array.Resize(ref Childs, Childs.Length+1); + Childs[Childs.Length-1] = child; + } + + //Recursively refreshes the current transform data + public void Refresh() { + //Local + if(Joint == null) { + Vector3 lp = Transform.localPosition; + Quaternion lr = Transform.localRotation; + LPX = lp.x; + LPY = lp.y; + LPZ = lp.z; + LRX = lr.x; + LRY = lr.y; + LRZ = lr.z; + LRW = lr.w; + } else { + XValue = Joint.X.GetTargetValue(true); + YValue = Joint.Y.GetTargetValue(true); + ZValue = Joint.Z.GetTargetValue(true); + Joint.ComputeLocalTransformation(XValue, YValue, ZValue, out LPX, out LPY, out LPZ, out LRX, out LRY, out LRZ, out LRW); + } + Vector3 ws = Transform.lossyScale; + WSX = ws.x; + WSY = ws.y; + WSZ = ws.z; + + //World + ComputeWorldTransformation(); + + //Feed Forward + foreach(Node child in Childs) { + child.Refresh(); + } + } + + //Updates local and world transform, and feeds the joint variable configuration forward to all childs + public void FeedForwardConfiguration(double[] configuration, bool updateWorld = false) { + //Assume no local update is required + bool updateLocal = false; + + if(XEnabled && configuration[XIndex] != XValue) { + XValue = configuration[XIndex]; + updateLocal = true; + } + if(YEnabled && configuration[YIndex] != YValue) { + YValue = configuration[YIndex]; + updateLocal = true; + } + if(ZEnabled && configuration[ZIndex] != ZValue) { + ZValue = configuration[ZIndex]; + updateLocal = true; + } + + //Only update local transformation if a joint value has changed + if(updateLocal) { + Joint.ComputeLocalTransformation(XValue, YValue, ZValue, out LPX, out LPY, out LPZ, out LRX, out LRY, out LRZ, out LRW); + updateWorld = true; + } + + //Only update world transformation if local transformation (in this or parent node) has changed + if(updateWorld) { + ComputeWorldTransformation(); + } + + //Feed forward the joint variable configuration + foreach(Node child in Childs) { + child.FeedForwardConfiguration(configuration, updateWorld); + } + } + + //Simulates a single transform modification while leaving the whole data structure unchanged + //Returns the resulting Cartesian posture transformations in the out values + public void SimulateModification( + double[] configuration + ) { + double[] px=Model.PX; double[] py=Model.PY; double[] pz=Model.PZ; + double[] rx=Model.RX; double[] ry=Model.RY; double[] rz=Model.RZ; double[] rw=Model.RW; + for(int i=0; i(); + if(segment != null) { + if(segment.Joint != null) { + segment.Joint.Remove(); + } + foreach(BioObjective objective in segment.Objectives) { + objective.Remove(); + } + Destroy(segment); + } + */ + + foreach(BioJoint joint in t.GetComponents()) { + joint.Erase(); + } + foreach(BioObjective objective in t.GetComponents()) { + objective.Erase(); + } + foreach(BioSegment segment in t.GetComponents()) { + Destroy(segment); + } + + for(int i=0; i stoppingDistance) { + //Accelerate + CurrentAcceleration = System.Math.Sign(CurrentError)*System.Math.Min(System.Math.Abs(CurrentError) / Time.deltaTime, maxAcc*Speedup); + } else { + //Deccelerate + if(CurrentError == 0.0) { + CurrentAcceleration = -System.Math.Sign(CurrentVelocity)* + System.Math.Min(System.Math.Abs(CurrentVelocity) / Time.deltaTime, maxAcc); + + } else { + CurrentAcceleration = -System.Math.Sign(CurrentVelocity)* + System.Math.Min( + System.Math.Min(System.Math.Abs(CurrentVelocity) / Time.deltaTime, maxAcc), + System.Math.Abs((CurrentVelocity*CurrentVelocity)/(2.0*CurrentError)) + ); + } + } + + //Compute new velocity + CurrentVelocity += CurrentAcceleration*Time.deltaTime; + + //Clamp velocity + if(CurrentVelocity > maxVel) { + CurrentVelocity = maxVel; + } + if(CurrentVelocity < -maxVel) { + CurrentVelocity = -maxVel; + } + + //Update Current Value + CurrentValue += CurrentVelocity*Time.deltaTime; + } + + public void SetEnabled(bool enabled) { + if(Enabled != enabled) { + Enabled = enabled; + Joint.Segment.Character.Refresh(); + } + } + + public bool IsEnabled() { + return Enabled; + } + + public void SetLowerLimit(double value) { + LowerLimit = System.Math.Min(0.0, value); + } + + public double GetLowerLimit(bool normalised = false) { + if(Constrained) { + if(normalised && Joint.JointType == JointType.Rotational) { + return Utility.Deg2Rad * LowerLimit; + } else { + return LowerLimit; + } + } else { + return double.MinValue; + } + } + + public void SetUpperLimit(double value) { + UpperLimit = System.Math.Max(0.0, value); + } + + public double GetUpperLimit(bool normalised = false) { + if(Constrained) { + if(normalised && Joint.JointType == JointType.Rotational) { + return Utility.Deg2Rad * UpperLimit; + } else { + return UpperLimit; + } + } else { + return double.MaxValue; + } + } + + public void SetTargetValue(double value, bool normalised = false) { + if(normalised && Joint.JointType == JointType.Rotational) { + value *= Utility.Rad2Deg; + } + if(Constrained) { + if(TargetValue > UpperLimit) { + value = UpperLimit; + } + if(TargetValue < LowerLimit) { + value = LowerLimit; + } + } + TargetValue = value; + } + + public double GetTargetValue(bool normalised = false) { + if(normalised && Joint.JointType == JointType.Rotational) { + return Utility.Deg2Rad * TargetValue; + } else { + return TargetValue; + } + } + + public double GetCurrentValue() { + return CurrentValue; + } + + public double GetCurrentError() { + return CurrentError; + } + + public double GetCurrentVelocity() { + return CurrentVelocity; + } + + public double GetCurrentAcceleration() { + return CurrentAcceleration; + } + + } + } + +} \ No newline at end of file diff --git a/Assets/External/BioIK/Setup/BioJoint.cs.meta b/Assets/External/BioIK/Setup/BioJoint.cs.meta new file mode 100644 index 00000000..9f8de9ce --- /dev/null +++ b/Assets/External/BioIK/Setup/BioJoint.cs.meta @@ -0,0 +1,19 @@ +fileFormatVersion: 2 +guid: 3531dc447309b4dfe8f7f4652bffd3fa +timeCreated: 1502289032 +licenseType: Store +MonoImporter: + serializedVersion: 2 + defaultReferences: [] + executionOrder: 0 + icon: {instanceID: 0} + userData: + assetBundleName: + assetBundleVariant: +AssetOrigin: + serializedVersion: 1 + productId: 67819 + packageName: Bio IK + packageVersion: 2.0d + assetPath: Assets/BioIK/Setup/BioJoint.cs + uploadId: 199125 diff --git a/Assets/External/BioIK/Setup/BioObjective.cs b/Assets/External/BioIK/Setup/BioObjective.cs new file mode 100644 index 00000000..40e3dace --- /dev/null +++ b/Assets/External/BioIK/Setup/BioObjective.cs @@ -0,0 +1,82 @@ +using UnityEngine; + +namespace BioIK { + + [AddComponentMenu("")] + public abstract class BioObjective : MonoBehaviour { + public BioSegment Segment; + public double Weight = 1.0; + + void Awake() { + + } + + void Start() { + + } + + void OnEnable() { + if(Segment != null) { + Segment.Character.Refresh(); + } + } + + void OnDisable() { + if(Segment != null) { + Segment.Character.Refresh(); + } + } + + void OnDestroy() { + + } + + public BioObjective Create(BioSegment segment) { + Segment = segment; + hideFlags = HideFlags.HideInInspector; + return this; + } + + public void Remove() { + for(int i=0; i 0) { + System.Array.Resize(ref Points, Points.Length-1); + } + } + } + + [System.Serializable] + public class DistancePoint { + public Transform Target; + public double Radius; + public double TPX, TPY, TPZ; + + public void SetTargetTransform(Transform t) { + Target = t; + if(Target != null) { + SetTargetPoint(Target.position); + } + } + + public void SetTargetPoint(Vector3 point) { + TPX = point.x; + TPY = point.y; + TPZ = point.z; + } + + public Vector3 GetTargetPoint() { + return new Vector3((float)TPX, (float)TPY, (float)TPZ); + } + + public void SetRadius(double radius) { + Radius = radius; + } + + public double GetRadius() { + return Radius; + } + } + +} \ No newline at end of file diff --git a/Assets/External/BioIK/Setup/Objectives/Distance.cs.meta b/Assets/External/BioIK/Setup/Objectives/Distance.cs.meta new file mode 100644 index 00000000..50b50180 --- /dev/null +++ b/Assets/External/BioIK/Setup/Objectives/Distance.cs.meta @@ -0,0 +1,19 @@ +fileFormatVersion: 2 +guid: 11a19f1e4ed2541b794d518fc3b60cce +timeCreated: 1502485934 +licenseType: Store +MonoImporter: + serializedVersion: 2 + defaultReferences: [] + executionOrder: 0 + icon: {instanceID: 0} + userData: + assetBundleName: + assetBundleVariant: +AssetOrigin: + serializedVersion: 1 + productId: 67819 + packageName: Bio IK + packageVersion: 2.0d + assetPath: Assets/BioIK/Setup/Objectives/Distance.cs + uploadId: 199125 diff --git a/Assets/External/BioIK/Setup/Objectives/JointValue.cs b/Assets/External/BioIK/Setup/Objectives/JointValue.cs new file mode 100644 index 00000000..a0688380 --- /dev/null +++ b/Assets/External/BioIK/Setup/Objectives/JointValue.cs @@ -0,0 +1,131 @@ +using UnityEngine; + +//!!!!!! +//This objective type is still under development +//!!!!!! + +namespace BioIK { + + //This objective aims to keep particular joint values acting as soft-constraints. Using this requires + //a joint added to the segment, and introduces some sort of stiffness controlled by the weight to the joint. + //Currently, you will require one objective for each >enabled< motion axis you wish to control. + [AddComponentMenu("")] + public class JointValue : BioObjective { + + public double TargetValue = 0.0; + public bool X, Y, Z; + + private BioJoint.Motion Motion; + private bool Valid; + private Model.MotionPtr Ptr; + private double NormalisedTargetValue; + + public override ObjectiveType GetObjectiveType() { + return ObjectiveType.JointValue; + } + + public override void UpdateData() { + if(Segment.Character.Evolution == null) { + return; + } + Valid = IsValid(); + if(!Valid) { + return; + } + if(X) { + Motion = Segment.Joint.X; + } + if(Y) { + Motion = Segment.Joint.Y; + } + if(Z) { + Motion = Segment.Joint.Z; + } + Ptr = Segment.Character.Evolution.GetModel().FindMotionPtr(Motion); + NormalisedTargetValue = Segment.Joint.JointType == JointType.Rotational ? NormalisedTargetValue = Utility.Deg2Rad * TargetValue : TargetValue; + } + + public override double ComputeLoss(double WPX, double WPY, double WPZ, double WRX, double WRY, double WRZ, double WRW, Model.Node node, double[] configuration) { + if(!Valid) { + return 0.0; + } + double loss = configuration[Ptr.Index] - NormalisedTargetValue; + return Weight * loss * loss; + } + + public override bool CheckConvergence(double WPX, double WPY, double WPZ, double WRX, double WRY, double WRZ, double WRW, Model.Node node, double[] configuration) { + return true; + } + + public override double ComputeValue(double WPX, double WPY, double WPZ, double WRX, double WRY, double WRZ, double WRW, Model.Node node, double[] configuration) { + if(!Valid) { + return 0.0; + } + return System.Math.Abs(configuration[Ptr.Index] - NormalisedTargetValue); + } + + public void SetTargetValue(double value) { + TargetValue = value; + } + + public double GetTargetValue() { + return TargetValue; + } + + public void SetXMotion(bool enabled) { + if(X != enabled) { + ResetMotion(); + X = enabled; + } + } + + public bool IsXMotion() { + return X; + } + + public void SetYMotion(bool enabled) { + if(Y != enabled) { + ResetMotion(); + Y = enabled; + } + } + + public bool IsYMotion() { + return Y; + } + + public void SetZMotion(bool enabled) { + if(Z != enabled) { + ResetMotion(); + Z = enabled; + } + } + + public bool IsZMotion() { + return Z; + } + + private void ResetMotion() { + X = false; Y = false; Z = false; + } + + private bool IsValid() { + if(Segment.Joint == null) { + return false; + } + if(!X && !Y && !Z) { + return false; + } + if(X && !Segment.Joint.X.IsEnabled()) { + return false; + } + if(Y && !Segment.Joint.Y.IsEnabled()) { + return false; + } + if(Z && !Segment.Joint.Z.IsEnabled()) { + return false; + } + return true; + } + } +} \ No newline at end of file diff --git a/Assets/External/BioIK/Setup/Objectives/JointValue.cs.meta b/Assets/External/BioIK/Setup/Objectives/JointValue.cs.meta new file mode 100644 index 00000000..b0f079b8 --- /dev/null +++ b/Assets/External/BioIK/Setup/Objectives/JointValue.cs.meta @@ -0,0 +1,19 @@ +fileFormatVersion: 2 +guid: dbc0a4a13428c4d64b8f01bf913c9c95 +timeCreated: 1502485934 +licenseType: Store +MonoImporter: + serializedVersion: 2 + defaultReferences: [] + executionOrder: 0 + icon: {instanceID: 0} + userData: + assetBundleName: + assetBundleVariant: +AssetOrigin: + serializedVersion: 1 + productId: 67819 + packageName: Bio IK + packageVersion: 2.0d + assetPath: Assets/BioIK/Setup/Objectives/JointValue.cs + uploadId: 199125 diff --git a/Assets/External/BioIK/Setup/Objectives/LookAt.cs b/Assets/External/BioIK/Setup/Objectives/LookAt.cs new file mode 100644 index 00000000..db9fefb0 --- /dev/null +++ b/Assets/External/BioIK/Setup/Objectives/LookAt.cs @@ -0,0 +1,123 @@ +using UnityEngine; + +namespace BioIK { + + //This objective aims to minimise the viewing distance between the transform and the target. + [AddComponentMenu("")] + public class LookAt : BioObjective { + + [SerializeField] private Transform Target; + [SerializeField] private double TPX, TPY, TPZ; + [SerializeField] private Vector3 ViewingDirection = Vector3.forward; + [SerializeField] private double MaximumError = 0.1; + + public override ObjectiveType GetObjectiveType() { + return ObjectiveType.LookAt; + } + + public override void UpdateData() { + if(Segment.Character.Evolution == null) { + return; + } + if(Target != null) { + Vector3 position = Target.position; + TPX = position.x; + TPY = position.y; + TPZ = position.z; + } + } + + public override double ComputeLoss(double WPX, double WPY, double WPZ, double WRX, double WRY, double WRZ, double WRW, Model.Node node, double[] configuration) { + double aX = 2.0 * ((0.5 - (WRY * WRY + WRZ * WRZ)) * ViewingDirection.x + (WRX * WRY - WRW * WRZ) * ViewingDirection.y + (WRX * WRZ + WRW * WRY) * ViewingDirection.z); + double aY = 2.0 * ((WRX * WRY + WRW * WRZ) * ViewingDirection.x + (0.5 - (WRX * WRX + WRZ * WRZ)) * ViewingDirection.y + (WRY * WRZ - WRW * WRX) * ViewingDirection.z); + double aZ = 2.0 * ((WRX * WRZ - WRW * WRY) * ViewingDirection.x + (WRY * WRZ + WRW * WRX) * ViewingDirection.y + (0.5 - (WRX * WRX + WRY * WRY)) * ViewingDirection.z); + double bX = TPX-WPX; + double bY = TPY-WPY; + double bZ = TPZ-WPZ; + double dot = aX*bX + aY*bY + aZ*bZ; + double len = System.Math.Sqrt(aX*aX + aY*aY + aZ*aZ) * System.Math.Sqrt(bX*bX + bY*bY + bZ*bZ); + double arg = dot/len; + if(arg > 1.0) { + arg = 1.0; + } else if(arg < -1.0) { + arg = -1.0; + } + double loss = System.Math.Acos(arg); + return Weight * loss * loss; + } + + public override bool CheckConvergence(double WPX, double WPY, double WPZ, double WRX, double WRY, double WRZ, double WRW, Model.Node node, double[] configuration) { + double aX = 2.0 * ((0.5 - (WRY * WRY + WRZ * WRZ)) * ViewingDirection.x + (WRX * WRY - WRW * WRZ) * ViewingDirection.y + (WRX * WRZ + WRW * WRY) * ViewingDirection.z); + double aY = 2.0 * ((WRX * WRY + WRW * WRZ) * ViewingDirection.x + (0.5 - (WRX * WRX + WRZ * WRZ)) * ViewingDirection.y + (WRY * WRZ - WRW * WRX) * ViewingDirection.z); + double aZ = 2.0 * ((WRX * WRZ - WRW * WRY) * ViewingDirection.x + (WRY * WRZ + WRW * WRX) * ViewingDirection.y + (0.5 - (WRX * WRX + WRY * WRY)) * ViewingDirection.z); + double bX = TPX-WPX; + double bY = TPY-WPY; + double bZ = TPZ-WPZ; + double dot = aX*bX + aY*bY + aZ*bZ; + double len = System.Math.Sqrt(aX*aX + aY*aY + aZ*aZ) * System.Math.Sqrt(bX*bX + bY*bY + bZ*bZ); + double arg = dot/len; + if(arg > 1.0) { + arg = 1.0; + } else if(arg < -1.0) { + arg = -1.0; + } + return System.Math.Acos(arg) <= Utility.Deg2Rad * MaximumError; + } + + public override double ComputeValue(double WPX, double WPY, double WPZ, double WRX, double WRY, double WRZ, double WRW, Model.Node node, double[] configuration) { + double aX = 2.0 * ((0.5 - (WRY * WRY + WRZ * WRZ)) * ViewingDirection.x + (WRX * WRY - WRW * WRZ) * ViewingDirection.y + (WRX * WRZ + WRW * WRY) * ViewingDirection.z); + double aY = 2.0 * ((WRX * WRY + WRW * WRZ) * ViewingDirection.x + (0.5 - (WRX * WRX + WRZ * WRZ)) * ViewingDirection.y + (WRY * WRZ - WRW * WRX) * ViewingDirection.z); + double aZ = 2.0 * ((WRX * WRZ - WRW * WRY) * ViewingDirection.x + (WRY * WRZ + WRW * WRX) * ViewingDirection.y + (0.5 - (WRX * WRX + WRY * WRY)) * ViewingDirection.z); + double bX = TPX-WPX; + double bY = TPY-WPY; + double bZ = TPZ-WPZ; + double dot = aX*bX + aY*bY + aZ*bZ; + double len = System.Math.Sqrt(aX*aX + aY*aY + aZ*aZ) * System.Math.Sqrt(bX*bX + bY*bY + bZ*bZ); + double arg = dot/len; + if(arg > 1.0) { + arg = 1.0; + } else if(arg < -1.0) { + arg = -1.0; + } + return Utility.Rad2Deg * System.Math.Acos(arg); + } + + public void SetTargetTransform(Transform target) { + Target = target; + if(target != null) { + SetTargetPosition(Target.position); + } + } + + public Transform GetTargetTransform() { + return Target; + } + + public void SetTargetPosition(Vector3 position) { + TPX = position.x; + TPY = position.y; + TPZ = position.z; + } + + public Vector3 GetTargetPosition() { + return new Vector3((float)TPX, (float)TPY, (float)TPZ); + } + + public void SetViewingDirection(Vector3 vector) { + ViewingDirection = vector; + } + + public Vector3 GetViewingDirection() { + return ViewingDirection; + } + + public void SetMaximumError(double degrees) { + MaximumError = degrees; + } + + public double GetMaximumError() { + return MaximumError; + } + } + +} \ No newline at end of file diff --git a/Assets/External/BioIK/Setup/Objectives/LookAt.cs.meta b/Assets/External/BioIK/Setup/Objectives/LookAt.cs.meta new file mode 100644 index 00000000..49a6c64e --- /dev/null +++ b/Assets/External/BioIK/Setup/Objectives/LookAt.cs.meta @@ -0,0 +1,19 @@ +fileFormatVersion: 2 +guid: f2576b795e2534cdbab732fca106cf79 +timeCreated: 1502416166 +licenseType: Store +MonoImporter: + serializedVersion: 2 + defaultReferences: [] + executionOrder: 0 + icon: {instanceID: 0} + userData: + assetBundleName: + assetBundleVariant: +AssetOrigin: + serializedVersion: 1 + productId: 67819 + packageName: Bio IK + packageVersion: 2.0d + assetPath: Assets/BioIK/Setup/Objectives/LookAt.cs + uploadId: 199125 diff --git a/Assets/External/BioIK/Setup/Objectives/Orientation.cs b/Assets/External/BioIK/Setup/Objectives/Orientation.cs new file mode 100644 index 00000000..b2f3fa5b --- /dev/null +++ b/Assets/External/BioIK/Setup/Objectives/Orientation.cs @@ -0,0 +1,106 @@ +using UnityEngine; + +namespace BioIK { + + //This objective aims to minimise the rotational distance between the transform and the target. + [AddComponentMenu("")] + public class Orientation : BioObjective { + + [SerializeField] private Transform Target; + [SerializeField] private double TRX, TRY, TRZ, TRW; + [SerializeField] private double MaximumError = 0.1; + + public override ObjectiveType GetObjectiveType() { + return ObjectiveType.Orientation; + } + + public override void UpdateData() { + if(Segment.Character.Evolution == null) { + return; + } + if(Target != null) { + Quaternion rotation = Target.rotation; + TRX = rotation.x; + TRY = rotation.y; + TRZ = rotation.z; + TRW = rotation.w; + } + } + + public override double ComputeLoss(double WPX, double WPY, double WPZ, double WRX, double WRY, double WRZ, double WRW, Model.Node node, double[] configuration) { + double d = WRX*TRX + WRY*TRY + WRZ*TRZ + WRW*TRW; + if(d < 0.0) { + d = -d; + if(d > 1.0) { + d = 1.0; + } + } else if(d > 1.0) { + d = 1.0; + } + double loss = 2.0 * System.Math.Acos(d); + return Weight * loss * loss; + } + + public override bool CheckConvergence(double WPX, double WPY, double WPZ, double WRX, double WRY, double WRZ, double WRW, Model.Node node, double[] configuration) { + double d = WRX*TRX + WRY*TRY + WRZ*TRZ + WRW*TRW; + if(d < 0.0) { + d = -d; + if(d > 1.0) { + d = 1.0; + } + } else if(d > 1.0) { + d = 1.0; + } + return 2.0 * System.Math.Acos(d) <= Utility.Deg2Rad * MaximumError; + } + + public override double ComputeValue(double WPX, double WPY, double WPZ, double WRX, double WRY, double WRZ, double WRW, Model.Node node, double[] configuration) { + double d = WRX*TRX + WRY*TRY + WRZ*TRZ + WRW*TRW; + if(d < 0.0) { + d = -d; + if(d > 1.0) { + d = 1.0; + } + } else if(d > 1.0) { + d = 1.0; + } + return Utility.Rad2Deg * 2.0 * System.Math.Acos(d); + } + + public void SetTargetTransform(Transform target) { + Target = target; + if(Target != null) { + SetTargetRotation(Target.rotation); + } + } + + public Transform GetTargetTransform() { + return Target; + } + + public void SetTargetRotation(Quaternion rotation) { + TRX = rotation.x; + TRY = rotation.y; + TRZ = rotation.z; + TRW = rotation.w; + } + + public void SetTargetRotation(Vector3 angles) { + SetTargetRotation(Quaternion.Euler(angles)); + } + + public Vector3 GetTargetRotattion() { + return new Quaternion((float)TRX, (float)TRY, (float)TRZ, (float)TRW).eulerAngles; + } + + public void SetMaximumError(double degrees) { + MaximumError = degrees; + } + + public double GetMaximumError() { + return MaximumError; + } + + } + +} \ No newline at end of file diff --git a/Assets/External/BioIK/Setup/Objectives/Orientation.cs.meta b/Assets/External/BioIK/Setup/Objectives/Orientation.cs.meta new file mode 100644 index 00000000..46333e20 --- /dev/null +++ b/Assets/External/BioIK/Setup/Objectives/Orientation.cs.meta @@ -0,0 +1,19 @@ +fileFormatVersion: 2 +guid: 819049a361f354ea4b9cac723854e0ff +timeCreated: 1502413797 +licenseType: Store +MonoImporter: + serializedVersion: 2 + defaultReferences: [] + executionOrder: 0 + icon: {instanceID: 0} + userData: + assetBundleName: + assetBundleVariant: +AssetOrigin: + serializedVersion: 1 + productId: 67819 + packageName: Bio IK + packageVersion: 2.0d + assetPath: Assets/BioIK/Setup/Objectives/Orientation.cs + uploadId: 199125 diff --git a/Assets/External/BioIK/Setup/Objectives/Position.cs b/Assets/External/BioIK/Setup/Objectives/Position.cs new file mode 100644 index 00000000..764de881 --- /dev/null +++ b/Assets/External/BioIK/Setup/Objectives/Position.cs @@ -0,0 +1,99 @@ +using UnityEngine; + +namespace BioIK { + + //This objective aims to minimise the translational distance between the transform and the target. + [AddComponentMenu("")] + public class Position : BioObjective { + + [SerializeField] private Transform Target; + [SerializeField] private double TPX, TPY, TPZ; + [SerializeField] private double MaximumError = 0.001; + + private double ChainLength; + private double Rescaling; + //private Vector3 Root; + + public override ObjectiveType GetObjectiveType() { + return ObjectiveType.Position; + } + + public override void UpdateData() { + if(Segment.Character.Evolution == null) { + return; + } + ChainLength = 0.0; + Transform[] chain = Segment.Character.Evolution.GetModel().FindObjectivePtr(this).Node.Chain;; + for(int i=0; i 0) { + SortByDistance(ref Hits); + for(int i=0; i 1.0) { + o = 1.0; + } + o = 2.0 * System.Math.Acos(o); + return Weight * 0.5 * (d*d + o*o); + */ + double d = Rescaling * ((TPX-WPX)*(TPX-WPX) + (TPY-WPY)*(TPY-WPY) + (TPZ-WPZ)*(TPZ-WPZ)); + double o = WRX*TRX + WRY*TRY + WRZ*TRZ + WRW*TRW; + if(o < 0.0) { + o = -o; + if(o > 1.0) { + o = 1.0; + } + } else if(o > 1.0) { + o = 1.0; + } + o = 2.0 * System.Math.Acos(o); + return Weight * 0.5 * (d*d + o*o); + } + + public override bool CheckConvergence(double WPX, double WPY, double WPZ, double WRX, double WRY, double WRZ, double WRW, Model.Node node, double[] configuration) { + double d = System.Math.Sqrt((TPX-WPX)*(TPX-WPX) + (TPY-WPY)*(TPY-WPY) + (TPZ-WPZ)*(TPZ-WPZ)); + double o = WRX*TRX + WRY*TRY + WRZ*TRZ + WRW*TRW; + if(o < 0.0) { + o = -o; + if(o > 1.0) { + o = 1.0; + } + } else if(o > 1.0) { + o = 1.0; + } + o = Utility.Rad2Deg * 2.0 * System.Math.Acos(o); + return d <= MaximumError && o <= MaximumError; + } + + public override double ComputeValue(double WPX, double WPY, double WPZ, double WRX, double WRY, double WRZ, double WRW, Model.Node node, double[] configuration) { + return 0.0; + } + + public void SetTargetTransform(Transform target) { + Target = target; + if(Target != null) { + SetTargetPosition(Target.position); + SetTargetRotation(Target.rotation); + } + } + + public Transform GetTargetTransform() { + return Target; + } + + public void SetTargetPosition(Vector3 position) { + TPX = position.x; + TPY = position.y; + TPZ = position.z; + } + + public Vector3 GetTargetPosition() { + return new Vector3((float)TPX, (float)TPY, (float)TPZ); + } + + public void SetTargetRotation(Quaternion rotation) { + TRX = rotation.x; + TRY = rotation.y; + TRZ = rotation.z; + TRW = rotation.w; + } + + public void SetTargetRotation(Vector3 angles) { + SetTargetRotation(Quaternion.Euler(angles)); + } + + public Vector3 GetTargetRotation() { + return new Quaternion((float)TRX, (float)TRY, (float)TRZ, (float)TRW).eulerAngles; + } + + public void SetMaximumError(double value) { + MaximumError = value; + } + + public double GetMaximumError() { + return MaximumError; + } + + public void SetNormal(Vector3 normal) { + Normal = normal; + } + + public Vector3 GetNormal() { + return Normal; + } + + public void SetSensitivity(float sensitivity) { + Sensitivity = Mathf.Clamp(sensitivity, 0f, 1f); + } + + public float GetSensitivity() { + return Sensitivity; + } + + public void SetLength(float length) { + Length = length; + } + + public float GetLength() { + return Length; + } + + } + +} \ No newline at end of file diff --git a/Assets/External/BioIK/Setup/Objectives/Projection.cs.meta b/Assets/External/BioIK/Setup/Objectives/Projection.cs.meta new file mode 100644 index 00000000..171e1019 --- /dev/null +++ b/Assets/External/BioIK/Setup/Objectives/Projection.cs.meta @@ -0,0 +1,19 @@ +fileFormatVersion: 2 +guid: fe8e654ecc18f42d3884dbbd9aed1624 +timeCreated: 1502485934 +licenseType: Store +MonoImporter: + serializedVersion: 2 + defaultReferences: [] + executionOrder: 0 + icon: {instanceID: 0} + userData: + assetBundleName: + assetBundleVariant: +AssetOrigin: + serializedVersion: 1 + productId: 67819 + packageName: Bio IK + packageVersion: 2.0d + assetPath: Assets/BioIK/Setup/Objectives/Projection.cs + uploadId: 199125 diff --git a/Assets/Scripts/HandTracking.meta b/Assets/Scripts/HandTracking.meta new file mode 100644 index 00000000..e3ee34f2 --- /dev/null +++ b/Assets/Scripts/HandTracking.meta @@ -0,0 +1,8 @@ +fileFormatVersion: 2 +guid: 2fd98f94acdd2294d93c612da351d805 +folderAsset: yes +DefaultImporter: + externalObjects: {} + userData: + assetBundleName: + assetBundleVariant: diff --git a/Assets/Scripts/HandTracking/HandTracking.meta b/Assets/Scripts/HandTracking/HandTracking.meta new file mode 100644 index 00000000..c5bffc22 --- /dev/null +++ b/Assets/Scripts/HandTracking/HandTracking.meta @@ -0,0 +1,8 @@ +fileFormatVersion: 2 +guid: 6535987dfca1ddd4c9d30a24447122eb +folderAsset: yes +DefaultImporter: + externalObjects: {} + userData: + assetBundleName: + assetBundleVariant: diff --git a/Assets/Scripts/HandTracking/HandTracking.unity b/Assets/Scripts/HandTracking/HandTracking.unity new file mode 100644 index 00000000..efe2ebcb --- /dev/null +++ b/Assets/Scripts/HandTracking/HandTracking.unity @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e338ca025d2bbdfc4af28b40904f2b38e461bea559978b4229fd1701529e8c69 +size 78811 diff --git a/Assets/Scripts/HandTracking/HandTracking.unity.meta b/Assets/Scripts/HandTracking/HandTracking.unity.meta new file mode 100644 index 00000000..2dd23e71 --- /dev/null +++ b/Assets/Scripts/HandTracking/HandTracking.unity.meta @@ -0,0 +1,7 @@ +fileFormatVersion: 2 +guid: 114fdda276c6fc147b9df81056c8a4b7 +DefaultImporter: + externalObjects: {} + userData: + assetBundleName: + assetBundleVariant: diff --git a/Assets/Scripts/HandTracking/HandTracking/GlobalVolumeProfile.asset b/Assets/Scripts/HandTracking/HandTracking/GlobalVolumeProfile.asset new file mode 100644 index 00000000..35fb6a3e --- /dev/null +++ b/Assets/Scripts/HandTracking/HandTracking/GlobalVolumeProfile.asset @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1cc2dbaf06bb8b671bb71a1ce3e0ebee43cf90575e84fb49cc2d63903112351b +size 2385 diff --git a/Assets/Scripts/HandTracking/HandTracking/GlobalVolumeProfile.asset.meta b/Assets/Scripts/HandTracking/HandTracking/GlobalVolumeProfile.asset.meta new file mode 100644 index 00000000..a5dad760 --- /dev/null +++ b/Assets/Scripts/HandTracking/HandTracking/GlobalVolumeProfile.asset.meta @@ -0,0 +1,8 @@ +fileFormatVersion: 2 +guid: 8ba096e0f19d4994ebcd28612e9c159d +NativeFormatImporter: + externalObjects: {} + mainObjectFileID: 11400000 + userData: + assetBundleName: + assetBundleVariant: diff --git a/Assets/Scripts/HandTracking/PuppetHand.fbx b/Assets/Scripts/HandTracking/PuppetHand.fbx new file mode 100644 index 00000000..7b38b90a --- /dev/null +++ b/Assets/Scripts/HandTracking/PuppetHand.fbx @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4c82b9fbaad1481d9dab80333b4aee1d6ddb56972ce8cb39dec5b3b82465ddff +size 382220 diff --git a/Assets/Scripts/HandTracking/PuppetHand.fbx.meta b/Assets/Scripts/HandTracking/PuppetHand.fbx.meta new file mode 100644 index 00000000..ec7f23bb --- /dev/null +++ b/Assets/Scripts/HandTracking/PuppetHand.fbx.meta @@ -0,0 +1,107 @@ +fileFormatVersion: 2 +guid: 3426073c46097c7408e65b9728d54018 +ModelImporter: + serializedVersion: 22200 + internalIDToNameTable: [] + externalObjects: {} + materials: + materialImportMode: 2 + materialName: 0 + materialSearch: 1 + materialLocation: 1 + animations: + legacyGenerateAnimations: 4 + bakeSimulation: 0 + resampleCurves: 1 + optimizeGameObjects: 0 + removeConstantScaleCurves: 0 + motionNodeName: + animationImportErrors: + animationImportWarnings: + animationRetargetingWarnings: + animationDoRetargetingWarnings: 0 + importAnimatedCustomProperties: 0 + importConstraints: 0 + animationCompression: 1 + animationRotationError: 0.5 + animationPositionError: 0.5 + animationScaleError: 0.5 + animationWrapMode: 0 + extraExposedTransformPaths: [] + extraUserProperties: [] + clipAnimations: [] + isReadable: 0 + meshes: + lODScreenPercentages: [] + globalScale: 1 + meshCompression: 0 + addColliders: 0 + useSRGBMaterialColor: 1 + sortHierarchyByName: 1 + importPhysicalCameras: 1 + importVisibility: 1 + importBlendShapes: 1 + importCameras: 1 + importLights: 1 + nodeNameCollisionStrategy: 1 + fileIdsGeneration: 2 + swapUVChannels: 0 + generateSecondaryUV: 0 + useFileUnits: 1 + keepQuads: 0 + weldVertices: 1 + bakeAxisConversion: 0 + preserveHierarchy: 0 + skinWeightsMode: 0 + maxBonesPerVertex: 4 + minBoneWeight: 0.001 + optimizeBones: 1 + meshOptimizationFlags: -1 + indexFormat: 0 + secondaryUVAngleDistortion: 8 + secondaryUVAreaDistortion: 15.000001 + secondaryUVHardAngle: 88 + secondaryUVMarginMethod: 1 + secondaryUVMinLightmapResolution: 40 + secondaryUVMinObjectScale: 1 + secondaryUVPackMargin: 4 + useFileScale: 1 + strictVertexDataChecks: 0 + tangentSpace: + normalSmoothAngle: 60 + normalImportMode: 0 + tangentImportMode: 3 + normalCalculationMode: 4 + legacyComputeAllNormalsFromSmoothingGroupsWhenMeshHasBlendShapes: 0 + blendShapeNormalImportMode: 1 + normalSmoothingSource: 0 + referencedClips: [] + importAnimation: 1 + humanDescription: + serializedVersion: 3 + human: [] + skeleton: [] + armTwist: 0.5 + foreArmTwist: 0.5 + upperLegTwist: 0.5 + legTwist: 0.5 + armStretch: 0.05 + legStretch: 0.05 + feetSpacing: 0 + globalScale: 1 + rootMotionBoneName: + hasTranslationDoF: 0 + hasExtraRoot: 0 + skeletonHasParents: 1 + lastHumanDescriptionAvatarSource: {instanceID: 0} + autoGenerateAvatarMappingIfUnspecified: 1 + animationType: 2 + humanoidOversampling: 1 + avatarSetup: 0 + addHumanoidExtraRootOnlyWhenUsingAvatar: 1 + importBlendShapeDeformPercent: 1 + remapMaterialsIfMaterialImportModeIsNone: 0 + additionalBone: 0 + userData: + assetBundleName: + assetBundleVariant: